NASA Contractor Report 3729 


Aeroelastic Analysis for Propellers 


Mathematical F ormulations and Program User’s M.anual 


R. L. Bielawa, S. A. Johnson, 

R. M. Chi, and S. T. Gangwani 

United Technologies Research Center 
East Hartford, Connecticut 


Prepared for 

Lewis Research Center 

under Contract NAS3-22753 


NASA 

National Aeronautics 
and Space Administration 

Scientific and Technical 
Information Branch 


1983 



PREFACE 


The G400PROP mathematical formulations and resulting computer code 
described herein were developed by United Technologies Research Center 
(UTRC) under contract NAS3-22753, "Development of a Comprehensive Aeroelastic 
Analysis for Propellers". This contract was through the Lewis Research Center 
of NASA with Mr. Oral Mehmed acting as contract monitor. The initial develop- 
ment of the G400 analysis was conducted at UTRC by Dr. Richard L. Bielawa 
under Corporate sponsored independent research and development. Extensive 
refinements to the analysis were made under sponsorship of the Langley Research 
Center of NASA and the U.S. Army Mobility R&D Laboratory, Langley Directorate 
as part of Contract NASI— 10960. Subsequently, further development was 
supported by Sikorsky and Hamilton Standard Divisions of United Technologies 
Corporation as well as the Structures Laboratory of the USRTL (AVRADCOM) under 
contract NAS1-16058. Sally Ann Johnson was responsible for developments in 
the expanded eigensolution, the PANPER interfacing and sample calculation 
portions of the report. She furthermore shared with Dr. Ray 11. Chi in the 
development of the unstalled unsteady airloads development. Dr. Santu T. 
Gangwani was responsible for the extentions of the UTRC stalled airloads 
theory to the higher subsonic Mach numbers . 
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Sally A. Johnson, 

Ray M. Chi, 

Santu T. Gangwani 

United Technologies Research Center 


SUMMARY 

Mathematical development is presented for a specialized propeller 
dedicated version of the United Technologies Corporation G400 Rotor Aero- 
elastic Analysis. This specialized analysis, G400PROP, simulates aeroelastic 
characteristics particular to propellers such as structural sweep, aerodynamic 
sweep and high subsonic unsteady airloads (both stalled and unstalled). 
Detailed formulations are presented for these expanded propeller related 
methodologies. Results are presented of limited application of the analysis 
to realistic blade configurations and operating conditions which include 
stable and unstable stall flutter test conditions. 


Sections are included for enhanced program user efficiency and expanded 
utilization. This material includes (1) a detailed description of the 
structuring of the G400PROP FORTRAN coding, (2) a detailed description of 
the required input data, (3) a detailed description of the output results, 
and (4) general information to facilitate operation and improve efficiency . 


* The research effort which led to the results in this report was financially 
supported by the NASA Lewis Research Center under contract no. NAS3-22753. 
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INTRODUCTION 


With recent renewed interest in propellers has come a redirection of 
design innovation and a need to advance the state-of-the-art in this 
technology. Major advances in materials and construction techniques on one 
hand, and the increasing ability to optimize basic aerodynamic efficiency 
on the other, are producing propeller designs which now require greater 
attention to analysis. A notable example is the so-called prop-fan which 
has extensive structural and aerodynamic sweep, relatively thin sections and 
operates in transonic flow conditions. This general acceleration in the growth 
of propeller state-of-the-art has especially increased the importance of 
structural integrity in a dynamic and/or aeroelastic environment. It is this 
aspect of propeller development to which the subject matter of this report is 
directed . 


Definition of Problem 


The assurance of satisfactory structural dynamic behavior and in particular 
aeroelastic (flutter) stability requires an accurate aeroelastic analysis 
specifically directed to the particular characteristics of advanced propeller 
designs. The purpose of this document is to describe the more important 
details of the G400PROP aeroelastic analysis developed to satisfy the analysis 
requirements of advanced propeller designs. 

The specific characteristics of advanced technology propellers as they 
relate to aeroelasticity can be readily identified. First, these propellers 
will generally continue to have sufficiently high aspect ratios, thereby 
justifying the treatment of them as beams for most aeroelastic problem areas. 

It is to be expected that, for some configurations with relatively low 
aspect ratio and high structural sweep, the beam theory formulation may have 
to be abandoned in favor of a more comprehensive plate theory formulation 
for some aeroelastic problem areas. The relatively low-cost advantages of 
beam theory formulations, however, together with the sustained applicability 
to configurations which do have high aspect ratios clearly justify develop- 
ment of a comprehensive aeroelastic analysis for propellers using beam theory. 

A second relevant characteristic of advanced technology propellers is 
the departure from the usual straight, torsionally rigid planforms. Within 
the technology available to build them, propellers are being designed with 
large sweeps and thinner sections to be rotated at significantly increased 
tip speeds in order to capitalize on the aerodynamic efficiencies which 
result. Structural sweep in a propeller blade is a relatively new and 
important aeroelastic consideration and clearly must be dealt with. The prime 
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importance of sweep is the large degree of coupling it introduces between 
bending and torsion. Based on experience with fixed wing sweep, this 
coupling must necessarily alter the aeroelastic behavior of propellers. 

The degree of aeroelastic involvement will also increase for these 
propellers due to the relative torsional softening caused by the thinner 
sections. Of particular significance is the increased susceptibility of 
these propellers to stall flutter, a condition usually experienced at the 
high thrust (high pitch angle) static flow conditions at take off. 

With the exception of structural sweep, the established validity of a 
beam formulation and the requirement to analyze known rotary wing aero- 
elastic phenomena (including stall flutter) form justifications for using 
a helicopter rotor aeroelastic analysis as a basis for advanced propeller 
aeroelastics. To this end, a highly successful helicopter and wind turbine 
aeroelastic analysis, the United Technologies Corporation (UTC) G400 Rotor 
Aeroelastic Analysis, was selected for enhanced development appropriate 
to advanced technology propellers. Two features of this analysis made it 
especially attractive to this application: First, it already had an 

advanced method for analyzing stall flutter, and second, it was formulated 
elastomechanically in a manner as to accommodate readily the inclusion of 
built-in structural sweep. Under NASA sponsorship, a development effort was 
therefore undertaken to modify a copy of the UTC G400 analysis into a compre- 
hensive aeroelastic analysis dedicated to the requirements of both general 
aviation and advanced technology propellers. 

Figure 1 presents an overview of this resulting computer code, G400PROP. 
The three principal types of inputs to the code are the physical description 
of the propeller (geometry, inertia, and elastic properties), the flight 
condition as defined by air density, speed of sound, propeller airspeed, and 
control (pitch) angle and, finally, an optional description of the detailed 
flow field resulting from nacelle blockage and wake/nacelle interaction 
considerations. The analysis generates dynamic equations of motion, which 
use a beam theory, normal modes basis and incorporates the higher order aero- 
elastic characteristics of structural sweep, structural twist, and unsteady 
airloads. For these dynamic equations, two principal solution types are 
produced: (1) eigensolutions , as defined in the Laplace Transform variable 

domain, and (2) time-history solutions appropriate to the calculation of 
transients. The principal uses of the eigensolution are the calculations 
of vacuum coupled modes (frequencies and mode shapes), and of those aero- 
elastic stability phenomena which can be readily linearized. The principal 
uses of the time-history solution are the calculations of transients result- 
ing either from strongly nonlinear aeroelastic stability phenomena or control 
inputs, of aerodynamic performance, and of harmonic responses of both hub 
loads and blade stresses resulting from harmonic aerodynamic excitation. 
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FLIGHT CONDITION 
(AIRSPEED. CONTROL ANGLE) 



Figure 1. Overview of G400PROP Aeroelastic Analysis 
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The version of the G400 used as a basis for the development of G400PROP 
was already structured in the form shown in Figure 1. However, several 
specific modifications and enhancements were required to develop G400PR0P. 
Principal areas in which modifications were required consisted of the following: 

1* All specific helicopter and wind turbine related modelings were 
to be stripped. 

2. The eigensolution portion of the analysis was to be upgraded 
to include both structural and aerodynamic sweep. 

3. The established UTRC unsteady stalled airloads methodology was 
to be extended to include high subsonic Mach number data. 

4. A differential equation (transfer function) modeling of unsteady 
unstalled airloads was to be developed and implemented. 

5. The analysis was to be configured to interact with an existing 
variable inflow analysis which describes the propeller /nacelle 
interactive flow field. 


Review of Existing Documentation 

The G400PROP propeller aeroelastic analysis described herein represents 
a specialization of an ongoing aeroelastic analysis development originally 
formulated for the unique aeroelastic characteristics of the composite 
bearingless rotor. It represented an advancement in the state-of-the-art 
with regard to the modeling of rotors with time-variable, nonlinear structural 
twist and multiple structural redundancy, as described in Reference 1. Since 
the publication of that report, the basic G400 program has evolved into a 
family of analyses with a completely general range of applicability in 
rotor type (articulated, hingeless, teetered) and vehicle application 
(helicopters, propellers, wind turbines). 

Most of the major documentation available on the G400 technology is 
contained in References 1 through 3. A review of their existing literature 
with regard to analytical alternatives to the G400 approach is contained 
in Reference 1 and a review of such alternate approaches now would be 
inappropriate. Reference 1 presents most of the basic ideas inherent in 
the G400 methodology upon which the present continuing development was made. 

The primary purpose of the present report, therefore, is to present the 
additional formulations required to meet the special requirements of analyzing 
advanced technology propellers. In addition, the development of this propeller 
oriented code coincidentally resulted in enhanced general program efficiency 
and capability. Therefore, a secondary purpose of this report is to increase 
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program user efficiency and to facilitate integration of G400PROP into 
industrial propeller design processes. 

Summary of New Technology 

All of the objectives identified for this development, as itemized in 
the first subsection above, were successfully met. The first major section 
to follow summarizes the unified theory developed for extending the original 
G400 twist related transformations to the more general case including 
structural sweep. The next four sections deal with elements of advanced 
aerodynamic modeling especially important to propellers. The first and fourth 
of these, which deal with aerodynamic sweep and inflow, respectively, are 
related more to aerodynamic ’'geometry" and are applicable even to steady 
flow conditions. The second and third of these aerodynamic sections deal 
with truly unsteady airloads and are those most directly related to propeller 
aeroelasticity . The following section presents the salient features of the 
enhanced eigensolution which include details of the modeling of perturbational 
sweep and unstalled unsteady related airloads. The following section presents 
details of limited application of the analysis to realistic propeller designs 
and appropriate operating conditions. The remaining four sections provide 
detailed program user information for the actual G400PROP computer program 
which implements the formulations presented herein. 
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LIST OF SYMBOLS AND FORTRAN EQUIVALENTS 


Symbol 

FORTRAN Equiv. 

Description 

A 

ARATE, AUNST 

Nondioensional aerodynamic section angle-of-attack time derivative 
(-ic/2U n ) 

[A] 

AA 

Inertia coupling matrix 

A 1* A 2* A 3 

CLA1 , CLA2 , CLA3 , 
(CMA1 y CMA2 , CMA3) 

Coefficients for inflow angle Pade augmented state variables, for either 
lift (or pitching moment) coefficients, (ND) 

a ,b ? c,d,e, f 

(None) 

Coefficients defining Pade approximant functionality 

a oL 

AO 

Aerodynamic section static lift curve slope, deg* 1 

a 

SPSD 

Sonic velocity at frees tream conditions, fps 

B 

TL 

Tip loss factor 

[B] 

BB 

Damping matrix for blade eigensolution 

VW B 4 

CLB1,CLB2, ...» 
(CMB1 , CMB2 , . . . ) 

Coefficients for pitching angle Pade augmented state variables, for either 
lift (or pitching moment) coefficients, (ND) 

b 

BL 

Number of blades 

[CJ 

CC 

Stiffness matrix for blade eigensolution 

c ,c ,c 

A w a 

(CADM.CARE), 
(CWDM, CWRE) , 
CATE 

Empirical coefficients multiplying A, and a, respectively, in dynamic 

stalled airloads functionality, (ND) 

C Ds’ C Ls’ C Ms 

CDS,CLS,CMS 

Aerodynamic section static drag, lift and pitching moment coefficients, 
respectively, as used in unsteady stalled airloads modeling 

c ,c ,c 

Du Lu Mu 

UNCD , UNCL , UNCM 

Aerodynamic section unsteady drag, lift and pitching moment coefficients, 
respectively, (ND) 

-C ,1C 
LI’ L2 

DELC1 ,DELC2 

Incremental lift coefficients (ND) 

c 

CHORD, CHORDB 

Blade section chord, ft and (ND) 

c d 

CD,ACD,CDTOT 

Section aerodynamic drag coefficient, (ND) 

c 0p’ c d» 

CDP , (CDS , CDSKNF) 

Section aerodynamic pressure and skin friction drag coefficients, respec- 
tively, (ND) 

c„ 

K 

CL, ACL 

Section aerodynamic lift coefficient, (ND) 

C V Cffl a 

(None) 

Unsteady aerodynamic section lift and pitching moment coefficients, respec- 
tively, due to pitching motion 

C ih* C «h 

(None) 

Unsteady aerodynamic section lift and pitching moment coefficients, respec- 
tively, due to plunging motion 

c 1 0 * c m 0 

(None) 

Steady-state section lift and pitching moment coefficients, respectively, to 
be used in conjunction with the incremental Pade section lift and moment 
coefficients, (ND) 

c m* c m c /4 

CM, ACM 

Section aerodynamic pitching moment about the quarter chord, (ND) 
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LIST OF SYMBOLS AND FORTRAN EOUIVS. (Cont'd) 


Symbol 

FORTRAN Equiv. 

Description 

c L 0 . c L r c L 2 

CLOPAN , CL1PAN , 
CL2PAN 

Quadratic components of the lift coefficient required by PANPER routine 

c D 0 . c Di.«D 2 

CDOPAN , CD1PAN , 
CD2PAN 

Quadratic components of the drag coefficient required by PANPER routine 

c u 

CRD PAN 

Blade section chord in the 4-coordinate system, (ND) 

Ac 

(None) 

Generalized incremental lift or pitching moment coefficient, defined by 
Pade theory, (ND) 


DCL.DCM 

Incremental section lift, or pitching moment, coefficients defined by 
Pade theory, (ND) 


(None) 

Section aerodynamic pressure and skin friction drags, respectively, lb/in. 

DUEAE^ , DUEAF ^ 

DUEAE, DUEAF 

Radial foreshortening of blade element point due to linear variations of 
k'th edgewise and i f th flatwise bending modes, respectively, (ND) 

DUEAO 

DUEAO 

Radial foreshortening of blade element point due to built-in sweep 

[DFDZ] 

DFDZ 

Partial derivative matrix of aerodynamic loadings with respect to inter- 
mediary perturbation vector 

EBi 

EB1B 

Nonlinear torsion stiffness parameter (to be multiplied by twist rate 
cubed), lb-ft^ 

eb 2 

EB2B 

Torsion to edgewise elastic coupling stiffness (to be multiplied by twist 
rate), lb-ft^ 

Ely, EI Z 

EIYB ,EIZB 

Section bending stiffness in flatwise and edgewise directions, respectivel> 
lb-in. ^ or (ND) 

e 

ER 

X 2 coordinate of coincident flat-lag hinge or hingeless blade offset 
point, in. 

IE 

(None) 

1st blade edgewise mode 

1F,2F,3F,4F,5F 

(None) 

1st through 5th blade flatwise modes 

GJ 

GJ,GJT,GJEFF 

2 

St. Venant (linear) torsion stiffness, lb-in. 

X E 

(None) 

4 

Edgewise area moment of inertia distribution, in. 

*F 

(None) 

4 

Flatwise area moment of inertia distribution, in. 

TV 

(None) 

Identify matrix of dimension m 

i 

(None) 

Square root of -1* 

k 

(None) 

Aerodynamic section reduced frequency (* uc/2U^) 

k A 

PRGS 

Area radius of gyration of tension carrying portion of blade section, in. 

k yio’ kz io 

KY10,KZ10 

Mass radii of gyration of blade section about axes through and perpendicular 
to the spanwise (X5) axis and in the chordwise and thicknesswise directions. 


respectively, in. 
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LIST OF SYMBOLS AND FORTRAN EOUIVS . (Cont'd) 


Symbol 

FORTRAN Equiv. 

Description 

L 

(None) , ELL 

Alternately (section aerodynamic lift distribution), and local radius of 
intersection of blade section midchord with section boundary, measured fror. 
axis <L , (ND) 

M 

MACH 

Mach number 

Mx 

(None) 

Section aerodynamic pitching moment distribution, in. -lb/in. 

m 

XMASSB 

Blade mass distribution (ND) 

m o 

XMO 

Reference blade mass distribution, taken to be that of the 5th blade seRir.eni 
lb-sec 2 /ft 2 

NEM,NFM,NTM 

(Same) 

Numbers of assumed flatwise, edgewise, and torsion natural M uncoupled M 
primitive modes, respectively 

n 

N 

Blade segment index 

n l »^2 * n 3 

(None) 

Unit vectors defining aerodynamically swept coordinate system, (ND) 

0( ) 

(None) 

Denotes order of magnitude 

P 

(None) 

Per rotor revolution, (ND) 

P 1 .P 2 .p 3 

P1,P2,P3 

Empirical coefficients for C Lu functionality 

P 

(None) 

Aerodynamically nondimensional Laplace transform variable 

p X5» p y 5 * p z 5 

SX5,SY5,SZ5 

Section shear load distributions in directions of axes in the 5-coordinate 
system, (ND) 

p l ,p 2 

CLP 1 , CLP 2 
(CMP1.CMP2) 

Yade poles for pitching, for either lift or pitching moment coef firl«-r*c 
(ND) 

p l ,p 2 

CLP1C.CLP2C, 

(CMP1C,CMP2C) 

Pade poles for plunging for either lift or pitching moment coefficients, 
(HD) 

Q lt Q 2 ,...Q 7 

Q1 »Q2 , . . .Q7 

Empirical coefficients for C Lu functionality 

{Q> 

Q 

Vector of blade degrees-of-f reedom 

q 

Q 

General expression for a response variable deflection 

qv k 

QVL(K) 

Blade k'th edgewise modal response variable 

q «i 

QWL(I) 

Blade i*th flatwise modal response variable 

q x 5 » q y 5 . q z 5 

XM5 , YM5 , ZMS 

Section moment load distributions about axes in the 5-coordinate svstem, 
(ND) 

q ej 

QTL(J) 

Blade j'th torsion modal response variable 

R 

R 

Rotor radius, ft 

R n 

(None) 

Reynolds number 

R 1* R 2* * * * R 8 

ARE1,ARE2, . . .ARE8 Empirical coefficients for Cjj u functionality 

r 

X 

Blade spanwise coordinate, measured from offset, e, in x $ direction (ND) 
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LIST OF SYMBOLS AND FORTRAN EQUIVS. (Cont’d) 


Symbol FORTRAN Equlv. Description 


Ar i ,Ax i 

DX , QUAD 

i'th blade spanwise segment (arc) length, (ND) 

s 

SSS 

Aerodynamic time 

s m 

SSSM 

Nondimens ional time measured from instant of stall onset 

T 

TENSB.TENST 

Tension at an arbitrary spanwise station, (ND) 

[TAS] 

TAS 

Coordinate transformation matrix relating "5" and "6 M coordinate systems, 
due to structural sweep, (ND) 

[TAS (A) ] 

ATAS 

Coordinate transformation matrix relating "5 M and "8 M coordinate systems, 
due to aerodynamic sweep, (ND) 


TR45 

Coordinate transformation matrix relating ”4" and "5" coordinate systems, 
(ND) 

1T.2T 

(None) 

1st and 2nd blade torsion modes 

t 

T 

Time, sec 

t, 

dm 

(None) 

Time when dynamic stall first occurs 

U 

U 

Total aerodynamic section inflow velocity, (ND) 

U N 

UN 

Vector sum of U T and U p aerodynamic section, inflow velocities, (ND) 

1r i * b’ p 

UR,UT,UP 

Total aerodynamic section inflow velocities in M 6" coordinate system radial, 
tangential and upflow directions, respectively, (ND) 

u r 5 » u t 5 * U P 5 

UR5QC,UT5QC, 

UF5QC 

Aerodynamic section inflow velocities in ”5” coordinate system radial, tan- 
gential and upflow directions, respectively, excluding blade motion, (ND) 

UELSET^ ,UELSFT^j , 

CE^SE k :.LI-ASF in 

UELSET,UELSFT, 

U ELA !v L , U l LA S F 

Radial foreshortening of blade element due to nonlinear variations of edge- 

flatwise and torsion mooes, <ND; 

u e 

UE 

Inward radial (x 5 ) foreshortening of blade element point due to combination 
of built-in sweep and elastic deformation, (ND) 

V T 

VEL.VELBAR 

Trimmed rotor flight speed, kts and (ND) 

w * 

WDIS 

Blade mass distribution, lb/in. 

L\\LV 

(None) 

Deflection correction terms due to second order twist effects, (ND) 

”e *^e 

VE.V7E 

Elastic deflections in the edgewise and flatwise directions, respectively, 
(ND) 

v y6 e ,Vz 6 e 

(None) 

Components of aerodynamic section inflow velocities in "6" coordinate syster 
tangential and upflow directions, respectively, due to blade motions, (ND) 

Av,Av 

(None) 

Deflection correction functions due to first order twise effects, (ND) 

Av r , Av , 

• v Bi 

uV EAj 

DVB,DVE,DVEA 

Edgewise motion deflection correction functions, (ND) 

v im 

VO,VINDND,PANVO 

Induced velocity based on momentum balance considerations, fps or (ND) 

fiW B k ’ iW e kj ’^EAj 

DWB , DWE , DWEA 

Flatwise motion deflection correction functions, (ND) 
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LIST OF 

SYMBOLS AND FORTRAN EOUIVS. (Cont'd) 

Symbol 

FORTRAN Equiv. 

Description 

X ,Y 
n n 

XALFW,YALFW 

Time variables used to calculate using a recursion formula 

x,y,z,w 

XN,YN,ZN,WN, Pade augmented state variables, for either lift (or pitching moment) 

(XNM,YNM,ZNH,WNM) 

x 5 ,y 5* z 5 

XX5.YY5.ZZ5 

Components of the 5-coordinate system, defined to be rotating with the hub, 
but at the blade coned and lagged position, (ND) 

x 4.y 4 ^A 

(None) 

Components of the 4-coordinate system, defined by vectors taken locally at 
the cyllndrically-oriented segment boundaries, (ND) 

* 4 LEB ,y4 LEB, 

Z4 leb 

X4LEB.Y4LEB, 

Z4LEB 

Components of the 4-coordinate system, taken at the intersection of the seg- 
ment leading edge with the segment boundary (ND) 

X4 QC’ y4 QC ,Z4 QC 

X4QC , Y4QC , Z4QC 

Components of the 4-coordinate system, taken at the intersection of the seg- 
ment quarter chord with the segment center, (ND) 

x 8*ys» 2 8 

(None) 

Components of the 8-coordinate system, defined similar to 6-coordinate sys- 
tem, but displaced by bending and oriented by the aerodynamic sweep angle 

y 6 * z 6 

YY6.ZZ6 

* 

Displacements in the 6-coordinate system, defined locally normal to built-ir. 
elastic axis 

y io * 2 io 

(None) 

Chordwise and thlcknesswise position coordinate, respectively, of an arbi- 
trary point within a blade section, (ND) 

y 10EA’*10EA 

Y10EA,Z10EA 

Built-in offset distances of elastic axis from x 5 axis, in edgewise and 
flatwise directions, respectively, (ND) 

^IOea-^IOea 

Y10EAC , Z10EAC 

Changes per segment (arc) length of built-in elastic axis edgewise and flat- 
wise offsets, respectively 

y 10cG 

Y10CG 

Chordwise distance of blade section mass center forward from the elastic 


axis, (ND) 


M 

o 

o 

y 10 3c /4 Y10QC,Y103QC 

Chordwise distances of blade section quarter chord and three quarter chord 
locations, respectively, forward from the elastic axis, (ND'* 

(Z) 

(None) 

Vector of intermediary aerodynamic variables 

a 

AL , ALF , ALTAB 

Section angle-of-attack, deg and rad 

a Dm 

ALDM 

Aerodynamic section dynamic stall angle-of-attack, deg 

°E 

ALTAB 

Effective aerodynamic section angle-of-attack, including effects of unsteady 
decay parameter, deg 

a RE 

ALRTCH 

Aerodynamic section reattachment angle-of-attack, deg 

°SS 

ALFSS 

Aerodynamic section static stall angle-of-attack, deg 

°TE 

(None) 

Section angle-of-attack when vortex nears the trailing edge 

a 

V 

ALFW 

Aerodynamic section unsteady decay parameter, rad 

a 

o 

(None) 

Mean angle-of-attack for oscillating airfoil, deg 

a 

(None) 

Angular amplitude for oscillating airfoil, deg 

Aa^ , Ac>2 

DELA1,DELA2 

Aerodynamic section angle-of-attack shifts to account for unsteady effects 


n 



v*> 


Symbol 


e 


6 B 




r ye * r 2e 

j j 

r * 

V 2e j 


Y 


Tv k 


' w i 


Y. 

C j 


6 < ) 


61*62 163*64 


n 


LIST OF SYMBOLS AND FORTRAN EOUIVS . (Cont'd) 

FORTRAN Eaulv . Description 

BETAT , FACPGT Alternatel y* to tal cone angle, rad, and Prandtl-Glauert transformation 

factor, (-/l-M 2 ) 


PREC,BETAB, 

BTAERO 

BETA1 

SHAPE 

GTNLY , GTNLZ 
GTNLYP , GTNLZP 

GAMMA 

GV 

GW 

GT 

(None) 

DELTAB , DELTAT 


Built-in Blade precone, deg or radians, as appropriate 

Empirical constant, normally equals 0.18 
Galerkin method integration weighting matrix 

Nonlinear j * th torsion modal weighting function for torsion excitation due 
to edgewise and flatwise force loadings, respectively, (ND) 

Nonlinear j f th torsion modal weighting functions for torsion excitation due 
to flatwise and edgewise moment loadings, respectively, (ND) 

Total aerodynamic sweep angle consisting of built-in and radial flow 
contributions, rad 

Deflection mode shape for the k' th edgewise normal mode, (ND) 

Deflection mode shape for the i ' th f latwise normal mode, (ND) 

Deflection mode shape for the j ' th torsion normal mode, (ND) 

Effective torsion mode shape due to integration of cosine components, (NL) 
Total blade lead angle, radians 


Denotes perturbational quantity 

PREL,DELTB , Built-in blade prelead, deg or rad, as appropriate 

DLAERO 

DELLC1 ,DELLC2 , Dynamic parameters used in dynamic stalled airloads functionality, (ND) 
DELLC3 ,DELLC4 


(None) The approximate distance measured from the G400 5-coordinate system segment 

boundary to the 4-coordinate system segment boundary, at a generalized blade 
chordwise location 


ETALE,ETATE 


n *n * • 
1 2 


ETA1 , ETA2 , . . . , 
ETA8 


The approximate distances measured from the G400 5-coordinate system segment 
boundary to the 4-coordinate system segment boundary, at the blade leading 
and trailing edges, respectively, (ND) 

Empirical coefficients for C Mu functionality 


6 


9 b 


e 

e 


e 


o 


6 


4 


TH 

THETA 

TW 

THE 

THAO, 

(TH75 ,THETA0) 
PCHPAN 


Total local blade pitch angle, radians 

Built-in blade pitch angle (structural twist), deg or rad 
Built-in twist rate, (ND) 

Elastic torsion deflection angle, radians 

Alternately, local static blade pitch angle, as defined by Pade filtering 
technique, or pitch angle due to input control angle, deg and rad 

Local blade pitch angle in the 4-coordinate system, radians 
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LIST OF SYMBOLS AND FORTRAN EQUIVS . (Cont'd) 

FORTRAN Equiv. Description 

THTIL Local blade perturbational pitch angle, as defined by Pade filtering tech- 

nique, rad 

CAPPA Aerodynamic pitch damping parameter, (ND) 

GME5 Structural sweep angle projection onto Xj-y^ plane, rad 

AGME5 Aerodynamic sweep angle projection onto X 5 «y 5 plane, rad 

GMF5 Structural sweep angle projection onto x^-z^ plane, rad 

AGMF5 Aerodynamic sweep angle projection onto x^-y^ plane, rad 

SWPLAM Aerodynamic section sweep angle of midchord measured in edgewise direction, 

(+) aft, deg 

A/. DGAM Difference in sweep angle between midchord and elastic axis, as measured 

in local chordwise direction, deg and rad 

X LAMBDA, (NONE) Alternately, aerodynamic rotor inflow, and eigenvalue (» o ♦ iui) 

(ST} XI Vector of excitations for the degrees-of-f reedom 

£ ZETA Angle which the radius measured from the hub to the intersection point of 

the section midchord with the section boundary, makes with the preconed 
and prelead-lagged feathering axis. 

{ ZETAMC Angle which the radius, measured from the hub to the intersection point of 

the section midchord with the section center, makes with the preconed and 
prelead-lagged feathering axis. 

RHO Air density, lb-sec^/ft^ 

SIGMA, ROOTE Alternately, rotor solidity, and real part of eigenvalue 

SIGXL,SIGYL, Integrations of Pade pole distributions for lift (and moment) 

SIGZL, SIGWL, 

(SIGXM, SIGYM, 

SIGZM.SIGWM) 

(NONE) Alternately, section thickness ratio, (• section thickness /semichord) , 

(ND) , and blade torsion stress, psi 

Total local blade inflow angle, radians 

Local blade static inflow angle, as defined by Pade filtering technique, 
deg or rad 

Local blade perturbational inflow angle, as defined by Pade' filtering 
technique, rad 

$ (NONE) Generalized Wagner function with compressibility corrections 

^ PSI (PSIREF) Blade azimuthal (angular) position, rad and (deg) 

£4 DP Nondlmensional time (azimuthal) step, rad 

u (NONE), ROOTI Alternately, frequency of motion, (Hz), and Imaginary part of eigenvalue 

Ww i* Wv k* W °J FF , EF , TF (Nondlmensional) uncoupled natural frequencies of i'th flatwise bending 

mode, k'th edgewise bending mode, and J f th torsion mode, respectively 


t PHI 

t PHIO 

o 

« PHTILN 
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LIST OF SYMBOLS AND FORTRAN EOUIVS. (Cont'd) 

Symbol FORTRAN Equlv. Description 

fl OMEGA, PRPM Rotor rotational frequency or speed (rpm) 


Subscripts 

( >a 

( )» 


( ) 

( >e 

< \ 

( ) 

( ) 

( ) 


EA 


LE 


LEB 


MC 


( ) 


MCB 


( > M 


( ) 


QC 


( ) 


TE 
( >u 

Superscripts 

() (A) 


( ) 


( ) 


00 


(M) 


(1) , v (2) 


( ) ,< ) 


( ) 

★ 

( ) 

( )' 

f ) 


Arising from aerodynamic loading 

Structurally built-in parameter, or conditions of blade immediately outboard 
of juncture 

Effects of dynamic origin 
Defined at the elastic axis 
Due to elastic deformation 
Pertaining to particular time step 
Defined at the blade leading edge 

Defined at the intersection of the section leading edge with the section 
boundary 

Defined at the intersection of the section midchord with the section 
center 

Defined at the intersection of the section midchord with the section 
boundary 

Components of vector quantities normal to segment midchord, or about 
local elastic axis 

Defined at the blade quarter chord 

Defined at the blade trailing edge 

Pertaining to unsteady stalled aerodynamic effects 


Pertaining to aerodynamic as opposed to structural 

Relating to Pade lift coefficient 

Relating to Pade pitching moment coefficient 

Pertain to first and second integrals defining the deflection correction 
function, respectively 

Nondimens ionalization by combinations of m 0 , R and/or ft 
Differentiation with respect to (Ot) 

Differentiation with respect to (r/R) 

Denotes evaluation at zero collective angle 
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STRUCTURAL TWIST AND SWEEP 


Principal Assumptions 

The aeroelastic analysis presented herein represents a specialization 
of the G400 analysis described in References 1 and 3. Since the publication 
of the initial documentation (Reference 1), various of the principal assump- 
tions enunciated therein have been either relaxed and/or extended. Most 
significant are the specific extensions which have been made to account for 
the more stringent modeling demands of high performance propellers (e.g. the 
prop-fan). Insofar as is practical within the scope of this report, the 
principal assumptions used herein are presented below: 

!• The rotor is rotating at a constant angular velocity, has infinite 
hub impedance, and is in steady translational flight. The orien- 
tation of the rotor in space is specified by appropriate Euler angles 
(pitch and roll). The orientation relative to the freestream is 
specified by means of a rotor angle-of-attack and a yaw angle. 

2. The elasto-meclianics of the blade are describeable within a beam 
theory framework with corrections of a kinematic nature to account 
for structural twist and sweep (see Figures 2a and b). The elements 
of beam theory analysis pertinent to the development, are the concepts 
of an elastic axis, the relationship of elastic bending to elastic 
torsion, and the definitions of these two elastic deformations. 

3. The elastic (torsion) axis is defined as the spanwise locus of 
shear centers of the two-dimensional blade (beam) sections taken 
perpendicular to this spanwise locus. Note that this definition 
treats the elastic axis as an abstracted section property, as 
contrasted with what one would measure in a bench test of an actual 
curved beam. In such a test, the locus of points where bending loads 
produce no torsion deflection (at the points of load application) 
would conform to the usual interpretation of the "elastic axis." 

This axis, however, would be different from the herein usage of 

the term to denote the locus of section shear centers. The built-in 
structural sweep (elastic axis offset), together with the elastic 
bending deflections, defines an elastic axis which is generally a 
space curve about which the local torsion deflections must take 
place. Thus, as shown in Figure 2a, each spanwise beam segment will 
not in general be defined parallel to the other segments. For the 
analysis of the beam-like elastic properties, the structurally 
swept blade (Figure 2a) is assumed to have its so-defined elastic 
axis "straightened out". This artificial straightening defines an 
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STRAIGHTENED ELASTIC AXIS 
WITH ARC LENGTH MAINTAINED 
CONSTANT 


b) EQUIVALENT BEAM FOR DEFINING BEAM ELASTOMECHANICS 

Figure 2. Basis for Use of Beam Theory for Structurally Swept Blade 




"equivalent beam" whose (straight) elastic axis has the same arc 
length as the original swept blade (see Figure 2b). 

4. The blade elasticity is described by the conventional (linear) beam 
bending and (nonlinear) bar torsion characteristics, as formulated 
by Houbolt and Brooks in Reference 4, for the above defined 
"equivalent beam". It is recognized that various deficiencies 
have been identified in these and other earlier formulations, 

both with respect to their adequacy for moderate to large bending 
deflections (References 5 and 6) and with respect to the proper 
modeling of pretwisted beams under tension (References 7-9)* 
However, there is not yet well established agreement either on the 
impact of these deficiencies on propeller elasticity or, more 
importantly, on a final proper reformulation. Thus, the continued 
use herein of the Houbolt and Brooks elastic formulations must be 
viewed as an eventually correctable deficiency of uncertain 
importance, to be addressed at some future date. 

5. The elastic bending and torsion deflections are "small" and 
respectively defined in a local sense normal to and along the space 
curve as defined by the built-in elastic axis. These deflections 
are defined as "small" in the sense that the elastic bending slopes 
and torsion deflection angles conform to the usual definition for 
"small" angles. 

6. The elastic bending and torsion deflections are describable using 
the "uncoupled" normal bending and torsion modes of the "equivalent 
beam". Thus, the deflections in the flatwise and edgewise direc- 
tions are respectively given by: 

NFM 

w e =Z y w(r)q w .(t) (la) 

i=i 1 1 


NEM 

v fl V" 


(lb) 


, and the elastic torsion deflections are given by: 

NTM 

0e=I Y e . (r)q 0 .(t) 

H J 1 


(lc) 
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Uncoupled modes are herein defined to be those beam vibration modes 
calculated assuming zero precone, prelag, pitch, twist, center-of- 
gravity offset and elastic axis offset. Thus y Wi > Y Vk > and . are 
all mutually uncoupled from each other. 

7. Blade elastic bending is defined by the conventional beam bending 
differential equations (as developed in Reference 4) wherein the 
usual independent spanwise variable is taken to be the arc length 
along the elastic axis. Furthermore, these bending differential 
equations are defined locally using the loadings locally normal 

to the built-in elastic axis. Within this context, explicit elastic 
bending-torsion coupling due to structural sweep is omitted in 
favor of implicit coupling due to inertial, aerodynamic and 
gravitational loadings taken with appropriate sweep related 
kinematics. 

8. The blade aerodynamic and structural twist distributions are non- 
linear. Additionally, the total (integrated) angle of structural 
twist is negligible beyond second order; cases of large local 
twist rates over short sections of span are not denied, however. 

9. Local radial foreshortening is defined relative to the equivalent 
beam defined in Figure 2. Contributions to radial foreshortening 
accrue from (a) the built-in structural sweep, i.e. that which 
restores the equivalent beam to the original swept planform (shown 
in Figure 2a) , (b) first order (linear) functions of bending, 
arising from built-in structural sweep, (c) second order (nonlinear) 
functions of bending each with elastic torsion arising from built-in 
structural sweep, and (d) second order functions each of both 
flatwise and edgewise bending. Note that this greatly relaxes 
principal assumption number 6 given in Reference 1. 

10. The elastic axis is coincident with the feathering (pitch) axis at 
the root of the blade. The built-in elastic is furthermore defined 
relative to the feathering axis. 

11. The blade flapping and lead-lag degrees-of -freedom used in Reference 1 
are assumed to be fixed at the built-in values. Thus, pinned root 
conditions are denied. It is to be noted that this assumption pertains 
only to the propeller dedicated aeroelastic analysis described herein. 

12. The blade distributions of center-of -gravity , aerodynamic center and 
center-of-tension (intersection of flatwise and edgewise neutral axes) 
are defined in two dimensional sections normal to the space curve 
elastic axis. These distributions are furthermore generally 
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noncoincident with the elastic axis and measured relative to the 
elastic axis. 

13. The blade sections have finite thicknesswise mass, but generally 
the thicknesswise location of the section center-of-gravity away 
from the chordwise principal axis is negligible. 

The above assumptions are used for the mechanical developments in the 
subsections which immediately follow. Assumptions regarding the basic imple- 
mentation of aerodynamic sweep are described in the next chapter. 


Basic Methodology for Structural Twist 

The present aeroelastic theory is characterized to a large extent by the 
kinematic modeling selected to describe the bending deflections of a pre- 
twisted, prebent beam (elastic axis taken as a space curve). The material 
presented draws heavily on the development of Reference 1 subject to above 
assumptions 5 through 10. 


As shown in Figure 3, the "5” coordinate system is defined by the preconed 
and prelead-lagged feathering axis. The "6" coordinate system is defined by 
unit vectors taken locally normal to the preswept elastic axis with the y 
direction arbitrarily taken parallel to the x - y plane. 6 


In the presence of only blade pitch angle, 6 , the "5" and "6" coordinate 
systems would be the same and the resulting "5" coordinate system deflections 
would be simple trigonometric resolutions of the flatwise and edgewise 
deflections as given by Equations (la) and (lb). With the addition of 
(arbitrary) structural twist, however, a simple trigonometric resolution 
transformation of flatwise and edgewise deflections is incapable of satisfying 
the beam "force" boundary conditions at the blade tip and an "integrated" 
trigonometric transformation is required. 


As is shown in Reference 1, the required integrated effect can be 
achieved by means of a trigonometric resolution transformation not on 
deflections, but instead on the second spanwise derivatives of the deflection, 
i.e. the curvatures : 

_//_// _ _// . ^ 

y 6 * v # cos®-w # sm® (2a) 


- » -// • ^ 
*6 : v e sin® 


+ we^cos® 


(2b) 
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This basic coordinate system transformation has the advantage that the 
aforementioned boundary conditions are always satisfied. In of themselves, 
however, trigonometric resolutions of curvatures are not directly useful for 
defining the blade kinematics. Instead, they serve as a starting point for 
deriving such a deflection based transformation. The required blade kinematic 
coordinate system transformation for deflections can then be derived from 
Equations (2) using integrations by parts and invocation of assumption 8. 

While the details of this integration procedure is straightforward, it is 
sufficiently tedious to be inappropriate to the scope of the present report. 
The integration yields a trigonometric resolution coordinate transformation 
in the usual form, as given in Reference 4 and elsewhere in the literature, 
but with the addition of various Reflection correction" functions due to 
twist : 


y 6 = (v e + Av- AV)cos@ - (w e - Aw- AW)sin® + 0(® /3 ) (3a) 

z 6 2 (v?e + Av- AV)sin® + (w e - Aw - AW)cos® + 0(© /3 ) (3b) 


where the underlined terms are, by assumption, negligible and where the 
deflection correction functions are defined by the following expressions: 


first order in twist: 


Av = dr + jf W w^df z dr l 

AW. jfVv.dP, 


(4a) 


(4b) 


second order in twist : 

AV* ^e/Awdr, ♦ Aw/ 2, dr 2 


(5a) 



(5b) 
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It is to be noted that the total twist rate, Q\ contains the built-in 
twist, 9g, and the time dependent elastic deflection, e; (- Ye • q 0 . ) . Thus , 
the deflection correction functions, Av, Aw, AV, and AW, nominally contain 
both linear and nonlinear combinations of the modal time variables (q w . , q^, 
and q Q ) . Specifically, the linear ones involve the built-in twist angle 
and are denoted with a "B" subscript. The nonlinear combinations involve 
the elastic twist angle and hence, are proportional to the product of 
with either q w ^ or q v ^. Herein, the products anc * are re ~ 

tained only in the Av and Aw (first order) correction functions and are 
denoted with an "e" subscript. The AV and AW (second order) correction 
terms retain only the contributions due to built-in twist rate, 0 fi , and 
hence, are strictly linear. Thus, the deflection correction functions de- 
fined generally in Equation (4) are given specifically as follows: 


f^irst_ order _in__twist : 

Av = Av b + Av t 

Aw = Aw b + Aw t 


where: 

HFU 

Av b = I Av- q w 

i«l 1 1 


NFM NTM 

Av,= I I Av e -q w ;q B . 

i.i j*i 'j 1 I 


NEM 


NEM NTM 

Aw,= I I An, q v q, 

M j»l I 1 


(6a) 

(6b) 

(7) 

( 8 ) 

(9) 

( 10 ) 
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second order in twist: 


AV = AV e (11a) 

AW = AW B (11b) 

where : 

NEM 

AV b =I AV b q (12) 

k = l * * 

NFM 

AW b =£ AW^q*, (13) 

i=l 


Additional deflection correction functions are defined with the consideration 
of structural sweep given in the following subsection. 


Kinematic Representation for Structural Sweep 
Approximations for Small Sweep 

Structural sweep is defined in a general sense wherein both chordwise 
and thickness offsets of the built-in elastic axis, y and z , 

EA 10 EA 

respectively, are admitted (see Figure 3). Within the context of the material 
in the preceeding subsection an approximation to the effects of structural 
sweep, for “small" values, can be obtained heuristically by considering the 
structural sweep to be “pre-bends” in the elastic axis. Within this context, 
the deflection correction functions defined by Equations (4) and (5) would 
be modified by the following substitutions: 


*w, q Wj ~ 

^ 2,0 EA 


X QWj — 

— Z 'k> CA 

(14a) 

w— 

-y.o EA 


r <\~ 

- *0 EA 
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where only those terms involving the elastic modal variables would be 
retained. For example, those terms involving only built-in twist and 
built-in sweep would be omitted. Thus, for small structural sweep 
Equations (6) and (11) would be modified as follows: 


A* — 

— A*+ Av ea 

Am — 

— Aw+ Aw tA 

m — 

— AV + AS7 ia 

AW — 

— AW + AW ca 


where: 



NTH 

I 

j" 





and where: 


= to d? . m 

AWE *i ■ I % * a: % 

O'. 


(14b) 


(15) 


(16) 


(17) 


(18) 


(19) 
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Coordinate Transformation Appropriate for Large Sweep 


For the general case of moderate to large structural sweep, the material 
in the preceeding subsection is inappropriate, but serves as an introduction 
to the material which follows. A successful modeling of structural sweep must 
account for the fact that the sweep angles defining elastic axis offset orien- 
tation are Euler angles and must be carefully defined. 


fif\d ^tj-of^plane_Positl : on \£ector_Comj>onent8 

The general modeling of the blade y and z kinematics due to combined 
structural twist and sweep is accomplished in the following steps: 


1. The elastic axis of the "equivalent beam" described in an above 

subsection is "distorted" back to the original planform defined by the 
built-in structural sweep and segment arc length distributions. 

This step essentially defines the position in space of the elastic 
axis space curve. This positioning requires the x , y and z offset 
distances of the centers of the segments as well as projections onto 
the x 5 ~y 5 and x,.-z,_ planes of the swept elastic axis line segments. 

These projections define the sweep angle distributions, A and A 
as shown in Figure 4. fi 5 ^5* 


2. The orientations of the elastic axis line segments define the local 
"6" coordinate system. x & is defined parallel to the axis of the 

elastic axis line segment; y is defined parallel to the x -y 

® 5 5 

plane, (+) in leading edge direction; z^ is orthogonal to x^ and yg, (+) 

in the normally positive thrusting motion. It should be stressed 
that the result of step 1 is to produce, in addition to the inplane 
and out— of— plane offsets (Ay^ and Az^) of the elastic axis from the 

(reference) x^ pitch axis, a radial foreshortening (Ax^) due to the 

constancy of the total arc length of the elastic axis. Mathematically, 
this Ax^ kinematic foreshortening is modeled differently and separately 

from the Ay^ and Az^ kinematic modeling. 


3. 


The blade segments of the blade configuration resulting from steps 1 
and 2 are then pitched and twisted about their respective elastic 
axis line segments (x & axis) to restore the blade back to its original 
built-in, but undeflected position. The pitch and twist angles for 

each segment are defined relative to the y axis. 

6 
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DISPLACED AIRFOIL 
SECTION DEFINING 
PITCH AND PRINCIPAL 
BENDING DIRECTIONS 



(FLATWISE) 
I W , Aw 

e 


(EDGEWISE) 


V ,Av 
e 


PITCH 

BEARING' 


FEATHERING 
(PITCH) AXIS 


ELASTIC AXIS 


BLADE SEGMENT 
ARC LENGTH, Ar 


(PARALLEL TO x & -y 5 PLANE) 


Figure 4. Euler Angles Defining Structural Sweep Transformation and Section Pitch 


LO 



4. The blade is then elastically deflected in torsion about the built-in 
space curve elastic axis to define a first set of "small" incremental 
y 5 an< ^ deflections. This first set of small incremental deflections 

is governed by Equations (15) through (19). 

5. The blade is then elastically deflected in flatwise and edgewise 
bending (in the presence of the torsion deflection) to define a 
second set of small incremental deflections. This second set of 
incremental deflections is measured in the "6" coordinate system 
and is governed by the basic deflection transformations defined by 
Equations (3) through (5). 


6. The second set of small incremental "6" coordinate system deflections 
defined in step 5 is transformed to the "5" coordinate system using 
a Euler angle transformation derived from sweep angle projections 
A and A , discussed in above step 1. 

C t- i- 


7. The results of steps 1, 4 and 6 are combined to define the and z 

position vector components. These results are summarized by the 
development which follows. 


First, the sweep angle projection distributions are defined using the 
built-in elastic axis line segment changes per segment length, the (invariant) 
segment arc lengths, A r , together with changes to the projection angles caused 
by elastic torsion deflection: 


A« 5 = sin 1 1 [( Av eaj ) - AV EA j 2) )cos® + (Aw ea <2) + AV%^ 2)/ )$in®j q^.| ( 20 ) 


A f 5 = sin ’' { + [“ ( + cos ® + ( av ea| 2>/ ) sin®] q*. } (2 1) 

wiiere Ay and A 2 ^ are the built-in changes per segment length. For 
EA EA ~ 

consistency with the definitions used for other previously defined radial 
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distributions these spanwise variable quantities are considered to be "derived" 
quantities calculated from the corresponding quantities defined in the 
chordwise and thicknesswise directions. Ay and Az , respectively. 

iU EA EA 

In practice, however, the "5" coordinate system quantities are the more 
accurately known and the "10" coordinate system quantities are derived using 
the "5" quantities according to: 


a Vk> ea = cos ^ + Az 5e A »"#. 


(22a) 


AZ '°EA = _ Ay 5EA Sin * B + A2 5 ea COS0 # 


(22b) 


Accordingly, Ay.„ and A z are input to the program and Ay and 
10 EA 10 EA EA 

Az are calculated internally using the inverse transformation of Equations 

5 ea 

( 22 ). 


The coordinate system transformation relating the pitch axis ("5") 
coordinate system with the swept ("6") coordinate system makes use of the 
sweep angle projections given in Equations (20) and (21): 


{*«} * M M 

{*»} • [tas''] M-MM 
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M * 


X 

- sin A es 

sin A e s 

X 

cosAf 5 

cos A* 

T s 

- X sin Afj 

sinAfsSinA^ 

cosAf 5 

cosAf s 


sinAf 


cosAf s 


(25) 


where : 


X-y/\- sin 2 Ae 5 - sin 2 Af 


(26) 


Equations (3), (15) through (19), (22), and (24) through (26) 
can then be combined to yield the required expressions for inplane and 
out-of-plane displacement: 


f y 5 1 f y '°EA C0S *B- Z I0 EA Sin0 B 1 

l z 5 J \ y io EA sin 0 B + Zio ea cos08 J 


f (Av ea .- AV ea .)cos 0 + (Aw ea . + AW EA .)sin© ^ 

j* * \^ v EAj~ ^EAjJ s i n ® “ (^ w EAj + AW EAj)COS® J 


! o 

(v # + Av - AV) COS© -(w*- Aw- AW) Sin 0 

(v # + Av - AV)sin® + (w # - Aw - AW)cos © 


(27) 


where : 


[«]■ 


0 I 
0 0 


:] 


(28) 
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and where v , w , Av, Aw, AV, AW are linear and nonlinear combinations of 
e e 

q , q and q , as per Equations (6) through (14). 

W. V. v 

i k j 

Radial_Pos it ion Vec_t_or_C£m£onent 

The modeling of the radial, x^, kinematics is accomplished in 

accordance with assumptions 3 and 9. The basis of the radial position 
modeling is the assumption of arc length constancy for each blade segment. 
Figure 5 presents a pictorial representation for a typical segment of the steps 
followed in this modeling process: 

1. The elastic axis of the "equivalent beam" segment is rotated to the 
built-in swept positions defined by the projection sweep angles, 

A and A , (given by Equations (20) and (21), with q *0). 

e S o 5 o e j 

The cosine foreshortening that results from this step is the 
built-in (constant) value and, for each segment, this "first" 
foreshortening is referred to as (dAx)^. 

2. The elastic axis orientation sweep projection angles are modified 
due to the elastic torsion deflections, in accordance with Equations 
(20) and (21). 

3. The results of step 2 (with the q dependency linearized) are 

combined with the elastic bending deflections, y^ and z^ , to 

e e 

produce the "second" foreshortening referred to as (dAx)^, as shown 

in Figure 5. This second foreshortening contribution is linear in 
the bending deflection variables, q and q but contains nonlinear 

W i V k 

(quadratic) combinations of these variables with the torsion 
deflection, q^ . 

4. Using the built-in projection sweep angles, A and A , the cosine 

®5 f 5 

o o 

foreshortening due to bending away from the undeflected elastic axis 
position is the "third" contribution to the foreshortening and is 
referred to as (dAx)^. This third foreshortening contribution is 
nonlinear in both the bending deflection variables, q and q . 

W i V k 
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Figure 5. Contributions to Incremental Radial Foreshortening Due to Sweep and 
Elastic Deformations 
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5. The results of steps 1 through 4 define the total incremental fore- 
shortening over an arbitrary blade segment, (dAx) . The total 
deflection for any segment is obtained by integrating the 
increments inboard of that location. These steps are summarized 
by the development given below: 

J5ui.lt^in_CorUr ^bution 

Reference to Figure 5 gives 


(dAx), = dr-dx = dr- ^/dr 2 - Ay| EA ~ Az 1 ea 


(29) 


= dr [i - y/~\ - ( Ay, 0EA /Ar) 2 - ( A z , 0ea / Ar ) 2 ] 


^°nt£ibu_t i£ n _Lin£ a £ in_Be.nd_ing_ 

The second foreshortening contribution is obtained by taking components 
of the bending deflection in the x__ direction: 


(dAx) 2 = dr[sinAf 5 cosA, 


A 


z 


/ 

6e “ sinAe 5 cosAf 5 



(30) 


where the ( ~ ) superscript denotes evaluation with zero collective angle since 
foreshortening relative to the pitch axis is invariant with collective angle. 
Each of the trigonometric functions of the sweep angle projections is then 
linearized with respect to q Q as are the bending shapes, 

j 

y and z . The details of this linearization are straightforward, but 
6 6 
e e 

sufficiently tedious to be beyond the intent of this report. After this 
operation, Equation (3) can be rewritten as: 


(dAx) 2 = [d(DUEAF, ) • q Wj + d(0UEAE H ) - q Vk 

+ d(UELSFTjj) • q*j Q®j + d(UELSET^.) • Qv k ^©j] 


(31) 
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Co n t£i^bu t_i£n_N£nl^irv ea r_in ji end inj> 

The third foreshortening contribution is obtained by calculating the 
cosine foreshortening and taking components in the direction: 


(dAx) 3 = co$A e5o CosAf 5o [i - ^/i- v e ' 2 - w« 2 ] dr 

s' cosA e5o cosAf^ ^ ( v4 2 + we ,2 )dr 


(32) 


< dAx >3 = cos Ae 5o cos Af 5o ' ^[ yv k y ^ q v k qv m + n q "« q "n] dr (33) 


where summations over i, k, m and n are implied, 

Combination of_Contr ibii t ions_ 

All of the three contributions to the incremental radial position 
deflection must be integrated: 


*5n = r n u «n 


(34) 


where r^ is the radial location of the nth segment of the "equivalent beam" 
and is the summation of segment arc-lengths up to the center of the nth segment. 
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and where 


u en = £ n [(dAx), + (dAx) 2 + (dAx) 3 ] (35) 

Symbolically, u is given by: 
e 

Ue * (DUEAO) + (DUEAF,)q w . 

+ (DUEAE k )- q Vk 

+ (UELSETk ) q v Qg + (UELSFT, )q w q a (36) 

t * J J 1 J 

+ j (UELASE km )q Vk q Vm + j (UELASF in ) q w .q Wn 
where variation in all terms with radial station is implied. 
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Torsion Excitation for Elastic Axis 
As A Space Curve 


As given in Reference 1, and as recognized elsewhere in the literature, 
the torsion differential equation is comprised of three basic parts. The 
first part consists of the usual elastic stiffening terms, and the second 
consists of combinations of distributed moment loadings. The third part is 
the wholly nonlinear torsion loadings accruing from distributed force 
loadings acting on moment arms provided by curvature in the elastic axis. 

As given in Reference 1, the torsion equation is given by: 


[gj 0e + + Ieb,!©' 2 - 6b)&'~ EB^e"] 




elastic stiffening 





V 



moment loadings 


(37) 


+ {*S\f!Hfr P 'S {r * )<ir 2~ /'P*5 (f 2 )dr 2 + < W r ' ) ] dr > 

I 1 

■ 2 »"jjfUlZ,' P « 5 <r * ldr2 " Z,' P »» (r2)dr 2‘ q '5 (r ' ) ] dr '}@ 


curvatures functions of force loadings 
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In Reference 1, the curvatures used in the (nonlinear) third portion 

of the torsion equation were assumed to arise entirely from the elastic bending 

deflections, and w”, as per Equation (2), As such, it can be shown 

that the nonlinear excitation term in Equation (37) can be reduced to the 

familiar difference of bending stiffness, AeI ( = El - EX ) , term: 

z y 



[(EI Z - EI y )ve w e "- (e A T + EB 2 ($b + ~ fle )0e) w e "] 


(38) 


This method for including the effect is attractive principally because 
of its simplicity and has been used to good advantage by numerous investiga- 
tors. Three difficulties exist with this method of implementation, however. 
The first difficulty relates to the fact that the implementation of Equation 
(38) is based on a "mode deflection" description of internal bending moment. 
The difficulty with a mode deflection formulation per se is two-fold. 

Studies of the characteristics of "mode deflection" (References 10 and 11) 
have established that convergence to accurate representations of internal 
bending moment is often not assured with a small number of modes. This 
accuracy problem is then compounded by the fact that the two components 
of this nonlinear excitation are subtractive. This is evidenced by the 
differencing of the section bending stiffnesses as indicated above. 

A second difficulty with using the AEI method relates to the assumed 
space curve character of the elastic axis. As such, torsion deflections 
are seen to contribute to inplane and out-of-plane deflections in the 
presence of bending (see Equations (8) and (10)). Thus, an analogous 
nonlinear excitation effect exists in both the flatwise and edgewise bending 
equations. In the framework of the G400 analysis, these nonlinear excita- 
tions in the bending equations are most practically implemented using a 
"force integration" approach. Consequently, the use of a AEI mode deflection 
implementation in the torsion equation together with a force integration 
implementation in the bending equations results in a (coupled) modal mass 
matrix which is generally nonsymmetric. A nonsymmetric mass matrix is not 
intrinsically a weakness for isolated rotor simulation and has been 
successfully used for years in that mode. However, the potential exists 
for spurious divergent response conditions caused by an inertia matrix 
becoming nonpositive-definite due to this deflection dependent nonsymmetry. 

The third difficulty with the Equation (38) formulation is that it is 
difficult to include the built-in curvature due to structural sweep. 

Equation (38) requires curvature information which is not generally available 
for the built-in geometry 
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Because of these difficulties, the conventional AEI approach of 
Equation (38) was abondoned in favor of a "force integration" approach. 
Accordingly, the Galerkin approach is applied to the nonlinear excitation 
term and integration by parts is used to achieve an intermediary step 
needed to eliminate the explicit curvature terms: 


/o'k r/; v= d ^ 


* p y 5 4 j£' y 0j Z 5 dr2<Jr i + (z s T + q y 5 ) 4 ^j y 5 dr ‘ 

" <y 5 /T - ^z 5 ) /o' ^ 2 5 dr,} dr 


(39) 


Since this term represents the nonlinear effects, it is reasonable to 
use a zeroth order approximation to the curvature terms wherein the structural 
sweep in assumed to be "small”. With this assumption, all the integrals in 
Equation (39) can be evaluated using the deflection correction functions 
defined in the above subsections. Thus, Equation (39) becomes: 


1 r *i dr r [ p y 5 C0S@ + P^ 00 ] 


" [ p z 5 cos@ " p y 5 sin@ ] 


(AO) 


+ [l(we~ Aw < 2, ~ Aw ,2) ) + q y5 cos®+ qj 5 sin®J 


~ ^yflj [ T ^ v e + Av (2) -Av (2 ’)- q^cos® + qy 5 sin®Jj’dr 
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where: 


r y flj 1 y 8j (We + z,o EA - Aw- AW) - ( Av ca .-AV ea .) (41a) 


T *Q\ y ^ v e + y,o EA + Av “ AV) “ ^ w EAj +■ AW ea .) ( 41b ) 



*,< we' ♦ z{o EA - Aw (2) '- AW* 21 ') - ( Av‘1 - Avi 2 /) 


J 


) — ( Av^ Aj - AV ea . 


(41c) 


?i g f Y e . (v' + y,' 0EA + Av (2)'- AV (2) ') - (Aw ea < 1 2) '+ AW ea J 2) '] 


(41d) 


Equation (40) represents the required form of the "force integration" 
implementation of the nonlinear torsion excitation term. To conclude this 
subsection, three observations can be made of the above formulation: 

1, Equations (41) all reduce to zero for zero structural sweep and 

zero elastic deflection, as would be expected from the behavior of 
Equation (38). 
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2. In Equation (40), the terms multiplying the nonlinear torsion 

weighting functions (r ye are actually the force and 

moment loadings defined for the linear excitations of the 
bending equations. The nonlinear torsion weighting functions. 
Equations (41), thus serve in effect, as the virtual deflection 
functions arising from torsion deflections appropriate to the 
bending generalized loads. 

3. The validity of the force integration approach is substantiated 

by the fact that the resulting terms in the torsion equation which 
represent rows of the inertia matrix (reflecting the integration 
of inertia forces) produce complete mass matrix symmetry and 
consequently insure positive-definiteness. 
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AERODYNAMIC SWEEP 


Principal Assumptions 


The unsteady airloads formulation incorporated in the G400PROP analysis 
is based on aerodynamic concepts originally developed for helicopter rotor 
blades. A characteristic of the aerodynamics of helicopter rotor blades 
is the generally large variability in local air velocities due to a 
combination of rotation with translational motion within the plane of the 
rotor. As a result, the aerodynamic formulations which have evolved are 
typically of a "strip theory" type with varying degrees of refinement to 
account for unsteady and swept flow effects. Such refinements typically 
are two-dimensional and applied in a heuristic manner based on the strip 
theory assumption. This is generally the approach followed herein. 

In addition to the basic strip theory assumption, the following specific 
sweep related assumptions are made. 

1. The local aerodynamic section sweep angle is defined by the angle 
the local airflow direction makes with the blade section taken normal 
to the midchord line (see Fig. 6). 

2. The section angle-of-attack is defined by the inflow and pitch 
angles measured within the section taken normal to the midchord 
line. 

3. For those cases wherein the "quasi-static" option is invoked, 

the effective angle-of-attack is defined (using above assumption 2) 
as the sum of the pitch and inflow angles. For this case, inflow 
angle is evaluated using local flow velocities at the 3/4 chord 
control point. 

4. For those cases wherein either of the specific, more advanced 
unsteady methods of the next two sections are invoked, the angle- 
of-attack or plunge variables are also defined using above 
assumption 2, but with inflow angle evaluation at the 1/4 chord 
control point. 

5. Airfoil drag is divided into two vectorial components (pressure drag, 
and skin friction drag) which are vectorially added to give the 
total drag. Pressure drag is that generally associated with compres- 
sibility and lift, and locally acts in the direction normal to 
midchord line, whereas skin-friction drag acts in the direction 

of the local flow velocity. 
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6. Skin friction drag varies with span (and hence Mach number) but 
is invariant with angle-of-attack. 

7. Lift, pitching moment, and pressure drag coefficients are determined 
by the angle-of -attack and Mach number measured in the section 
normal to the midchord line. The lift, pitching moment, and pressure 
drag are determined by their so-defined coefficients and the 
dynamic pressure based on the velocity components normal to the 
midchord line. 

8. Skin friction drag is determined by the Mach number and dynamic 
pressure based upon the total vector sum of all components of the 
local total flow. 

These assumptions form the basis of the development which follows. 


Basic Modeling Characteristics 
Aerodynamic Sweep Angle Projections 

As was developed in the previous section, the appropriate axis for 
defining structural sweep is the elastic axis. The appropriate axis for 
defining aerodynamic sweep, however, is, by assumption 1, the locus of 
midchords. Therefore, the aerodynamic sweep angles are defined as: 


A ( iJ = Ae 5 + (COS 8 AA-yg e ) 


(42) 


A* = A f + (- sin 8 AA+Zg # ) 

5 5 


(43) 


where AA is the difference in sweep angle between the midchord and the 
elastic axis, as measured in the local chordwise direction, and can reasonably 
be assumed to be a "small" angle. Since the elastic bending slopes of 
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y and z , are also assumed to be small, the small angle assumption can 
o o 

e e 

safely be used on the parenthetical terms in Equations (42) and (43). 


Similar to the coordinate system transformation developed in the 
previous section (Equations (23) through (26)) and pictorially defined in 
Figure 4, a coordinate system transformation can be formulated for aerodynamic 
sweep. In particular, such an aerodynamic sweep transformation is needed 
to relate velocity and loading components in the "5" coordinate system to 
those quantities in the n 8 M coordinate systems. The "8" coordinate system 
is defined similarly to the "6" coordinate system but is additionally 
displaced by the elastic bending deflections, and is rotated relative to 
the "5" coordinate system by the aerodynamic sweep angles, 

(A) (A) 

A and A , given by Equations (42) and (43), respectively. Using the 
e 5 5 

vector decompositions shown in Figure 7, the following definitions are made 
for an arbitrary spanwlse element: 


n^ = unit vector along the span of the midchord of the blade element 

n 2 = unit vector perpendicular to 3^ and parallel to the x -y^ plane, 
positive forward 

n^ = unit vector mutually perpendicular to both and If , righthand rule. 


Similar to Equation (24) the following relationship can be written: 



( 44 ) 


(A) 

where the elements of [TAS ] are defined using Equations (fc5) and (26), 
but with the aerodynamic sweep angles. 


Velocity Decompositions 

As shown in Figure 7, the local velocity vector including only environ- 
mental effects (i.e., neglecting blade motion) can be expressed in either 
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Figure 7. Pictorial Representations of Coordinate Systems and Velocity Decomposition 
Relating to Aerodynamic Sweep 
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coordinate system: 


u= Ur 5 T- Ut 5 T + Up 5 k 

= UR 0 it, - Uj^ + Up 0 n 3 


( 45 ) 


Use of the above coordinate system transformation together with the local 
velocities due to blade motion, v and v , yields the following useful 
fsrm: y 6 Z 6 



( 46 ) 


Note the U , U , and U can include not only the components due to 
R 5 T 5 P 5 

axial flow and variable inflow distortion, as discussed in a later section, 
but any perturbation which might accrue from hub motion. Consideration of 
such hub motion is presently outside the scope of the study, however. 

Other relationships which are needed to formulate the airloading distri- 
bution are the magnitude of the total velocity, U, and the component normal 

to the midchord, U : 

* N 


u = (sgn Ut) 



2 2 2 
Uj + Up + Ur 


( 47 ) 


u N = (sgn Ur) + Up 


( 48 ) 
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Airload Distributions 


Using standard strip theory techniques, the local lift, pitching moment, 
pressure drag, and skin friction drag distributions can, respectively, be 
written as: 


L * P-C- ArU* c x (aNtM N ) 


(49a) 


^2 c2 ArU(5 Cm c/4 (aN,Mn) 


(49b) 


Dp= ^2 CArU N c<j P(CIn» m n) 


(49c) 


D s = /5f cAr U 2c d s ( M ) (49d) 

where : 

a N = + tan“‘ (Up/Ur) (50) 


m n = u n /a® 


(51) 


M = U/0» 


(52) 


c dp s C d^ a N> M N^“ C «J s ^ M h) 


(53) 


The pitch angle seen in the "8" coordinate system, 0 N , is obtained from 
the nominal pitch angles with consideration of the integrated effect of the 
cosine components: 
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( 54 ) 


NTM 

9 N Z Q 0 + 9 0 + £ Ye Qft 

i= I J * 


where: 




2 [ x ®j + f 0 y «j (cos *%,> - i )dr,]cos A ( 
+ [jf X$. sin A <A) (r,)dr,Jsin A* a) 


and where: 


A< a > = sin 


-i 


Ar 


+ AA 


(55) 


(56) 


Comp£nents_in the ^8^ c^oo^rd^inajt e_S£S tern 

The airload distributions given by Equations (49) are then resolved 
to the "8" coordinate system using the components of the inflow velocity: 


p Ox 0 = ^jCAr UU R C,j s 


(57) 


p 0 y 8 s -^fCAr[u N (C dp Ur- C A Up) -►UUtC^] (58) 

P a Z8 = ^2 CAr [ U n( c £U T + C dp Up) + UUp Cds ] (59) 

q ax 8 s ^SC?Ar[uJ c mc/4 - u*e] + y, 0c/4 (p Oz# COS0 - FJ, sine) (60) 
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and where: 


using advanced unsteady airloads 
option 

(61) 


K= -j- c(sgn U) 


• -■F-yio 


C ,lo C/4 


forward flow, quasi-static 
option 


C ^ l0 C/4 


: reversed flow f quasi-static 
option 


£ompon£nts_in _the ^ojor^iM_te_s^sjtem 

Above Equations (57) through (60) define the airload components 
with directions aligned with the deflected blade segments. The appropriate 
airloads needed are those defined in the "6", or undeflected coordinate 
system: 

P % = P % + “ P a !e ( ' sin9 + z «> (62) 



■v 


P„ (COS0AA - yJ) 

°*8 6t 


(63) 


p » t6 = p »i, + p »«, ( ‘ 5ir ' 0A ' l+2 *« ) 


(64) 


q °* 6 : q °*8 


(65) 


48 



Remarks Concerning Application 


Above Equations (63) through (65) define the most general information 
of the airload distributions to be used in the bending equations and in the 
nonlinear excitation term of the torsion equation (see Equation (41)). 

These expressions for airloadings are very nonlinear and, in the above form, 
are only suitable for utilization in the time-history solution. For eigen- 
solution purposes, they must be completely expanded to yield all explicit 
linearized perturbations of the modal variables and their derivatives, 

* 

6q , 6q , 6q , 6q ... The details of this perturbational expansion are 

W V 0 w 

i k j i 

straightforward, but sufficiently tedious to warrant excluding them from the 
present report. 
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UNSTEADY STALLED AIRLOADS 


A detailed analysis of dynamic stall experiments has led to a semi- 
analytic methodology characterized by a set of relatively compact analytical 
expressions, called synthesized unsteady airfoil data, which accurately 
describe in the time-domain the unsteady aerodynamic characteristics of 
stalled airfoils (Reference 12). Under the present study, the unsteady 
stalled airloads methodology was expanded for propeller applications by 
synthesizing similar unsteady loops at subsonic Mach numbers which are 
higher, more relevant than those used in the earlier study. More specifically, 
the high Mach number data contained in References 13 and 14 were reduced to 
synthesized form within the established Reference 12 framework. 


Review of Basic Methodology 


Dynamic Stall Model 

The analytical model of dynamic stall, described herein, includes the 
main physical features of the dynamic stall phenomenon as observed in 
oscillationing airfoil tests. A brief description of dynamic stall events 
is given below. 

When an airfoil experiences an unsteady increase in angle-of-attack 
beyond the static stall angle, a vortex starts to grow near the leading edge 
region. As the angle continues to increase, the vortex detaches from the 
leading edge and is convected downstream near the surface. These events are 
shown schematically in Figure 8. The suction associated with the vortex 
normally causes an initial increase in lift. The magnitude of the increase 
depends on the strength of the vortex and its distance from the surface. 

The streamwise movement of the vortex depends on the airfoil shape and the 
pitch rate. The relative distance between the vortex and the airfoil varies 
according to the kinematics of the airfoil. That is, it depends on charac- 
teristics such as the pitch rate and the instantaneous angle-of-attack. 

As the vortex leaves the trailing edge, a peak negative pitching moment is 
obtained. The airfoil then remains stalled until the angle-of-attack drops 
sufficiently so that reattachment of the flow can occur. The present method 
incorporates all of these events. For example, the strength of the vortex 
is made a function of the angle when the vortex leaves the leading edge 
(moment stall angle). The higher the moment stall angle, the higher the 
strength of the vortex. 
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VORTEX FORMATION NEAR 
LEADING EDGE 



MOVEMENT OF THE VORTEX 
OVER THE AIRFOIL 


Figure 8. Dynamic Stall Modeling 
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Parameters Influencing Dynamic Stall 


The unsteady lift, drag, and pitching moment coefficients of the airfoils 
obtained from the two-dimensional oscillating airfoil tests show a large 
degree of hysteresis when plotted as functions of angle-of-attack, 
particularly when the reduced frequency and the maximum angle-of-attack are 
sufficiently high. Figures 9a, 9b, and 9c show an example of typical 
loop data obtained from the oscillating airfoil test. The amount of 
hysteresis and the shape of the loops vary in a highly nonlinear fashion with 
such test parameters as amplitude, mean angle, and reduced frequency. 

The results of the oscillation airfoil tests clearly indicate that the 

dynamic characteristics of an airfoil depend on the following main parameters: 

(1) airfoil shape and sweep; (2) Mach number; (3) Reynolds number; (4) reduced 

frequency, k; (5) oscillation amplitude, a ; and (6) mean angle-of-attack, a . 

o 

The first three of these parameters affect both the static and the 
dynamic characteristics of the airfoil, while the last three parameters 
represent purely dynamic parameters. Since most rotor aeroelastic analyses 
employ time-history solution techniques for computation of the aerodynamic 
loading acting on the rotor blades , frequency domain parameters such as reduced 
frequency or amplitude, etc., are inappropriate for use in these time domain 
simulations. Moreover, for arbitrary motion it is difficult to describe the 
reduced frequency, the amplitude of oscillation, or the mean angle-of-attack 
of a rotor blade section in a precise manner. As a result, an alternative 
set of dynamic parameters, which are appropriate for the time domain simula- 
tions, is defined. The parameters replacing k, a, and a in the present method 
are: (4) the instantaneous angle-of-attack, a; (5) the nondimen sional 

pitch rate, A; and (6) the unsteady decay parameter, a , which accounts for 
the time history effects of the change in a, and is based upon the Wagner 
function. 

For the sinusoidally oscillating airfoil, these three parameters can be 

feasily expressed in terms of the reduced frequency, the amplitude, and the 

mean angle-of-attack. Also, they can be easily evaluated for rotor blade 

sections in a stepwise manner and are very convenient to use for the prediction 

of the onset of dynamic stall and for the determination of the unsteady airloads. 

Thus, the present method determines, through the synthesization process, 

the effect of these selected parameters ( a. A, a ) on the dynamic stall 

w 

characteristics of the airfoils by utilizing the data from the oscillating 
airfoil tests. The synthesization process used herein essentially involves 
curve-fitting of the test loop data to the prescribed analytical expressions, 
with the objective of determining the unknown parameters or coefficients 
embedded in the analytical expressions. The analytical expressions are 
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UNSTEADY LIFT COEFFICIENT, C, 



ANGLE OF ATTACK, a (deg) 


Figure 9a. Typical Unsteady Lift Coefficient Loop Data, SC 1095 Airfoil, M =0.3, 
a Q = 12.0 deg, a = 8.0 deg, k = 0.1 
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UNSTEADY PITCHING MOMENT COEFFICIENT, C 





UNSTEADY DRAG COEFFICIENT, C 



Figure 9c. Typical Unsteady Drag Coefficient Loop Data, SC 1095 Airfoil, M =0.3, 
a Q = 12.0 deg, a = 8.0 deg, k = 0.1 


obtained mostly by mathematical or empirical means and in general they 
represent simple quantitative approximations to the various observed 
physical features of the dynamic stall phenomenon. 

Definition of the Unsteady Decay Parameter, a 
w 

For a two-dimensional airfoil going through an arbitrary change in 

angle-of-attack, one can describe an instantaneous effective angle-of-attack, 

a , by using Duhamel's integral (Reference 10) as given below: 

E 


a E (s) = a (0) 4> c (s,m) + J 4 s - /9 4> c a ‘» M) dcr 

ft * ^ 


( 66 ) 


where a(0) corresponds to the initial angle-of-attack, M represents Mach 

number, <j> (s,M) is the response to step change in a (a compressibility 
c 

corrected form of the Wagner function), and s is the nondimen sional time 
as given by: 

» * § X' Ud, ‘ <67) 


The unsteady decay parameter. 


ot , to be used extensively in the present 


method, is defined as follows: 


a w = a (s)-a E (s) 


( 68 ) 


The a parameter physically represents the difference between the instan- 
w 

taneous angle, a and the effective angle, a , and therefore accounts for the 

E 

time-history effects of the change in a. This physical description of a 

w 

is valid for attached flow conditions only. In the present method, the a 

parameter is most useful for predicting the onset of dynamic stall, and 
for convenience, it is also used to describe approximately the unsteady 
coefficients after the stall. 
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The effects of compressibility are incorporated in the definitions of 
<* w by the use of the general or compressible Wagner function (see also 

Reference 15) obtained from the following approximate relationship 


(69) 

4> c (s,M)« [l.o -O.I65e -00458 * (, " m2) -0.335 e _0,3s /y.-M* 

Computation of Dynamic Parameters 

For the sinusoidally oscillating airfoil, where the motion of the airfoil 
is completely known, the parameters a, A, and a can be obtained analyticaly 

as given below: 

a = 5 sin ks (70) , 

A = kS cos ks (71) 

a w = y, (k,M) k5 cos ks + y 2 (k,M) ff sin ks (72) 


where k, s, and M represent reduced frequency, nondimensional time, and 
Mach number, respectively. The y and y functions are described by: 


y (k M) - 0-165 (l-M z )(0.0455) , 0,335 (>-M )(0.3) 

' * k 2 +0-M 2 ) 2 (0.0455) 2 k 2 +(l-M 2 ) 2 (0.3) 2 


(73) 


r z (k, m) = 


0.1 65 k 2 

k 2 + (l-M 2 ) 2 (0.0455) 2 


+ 


0.335k 2 

k 2 + H-M 2 ) 2 (0.3) 2 


(74) 
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In contrast to the closed form evaluations obtainable for sinusoidal 
motion, numerical evaluations of these three section dynamic parameters must 
be obtained for arbitrary motion in the time domain. This is accomplished in 
a stepwise manner utilizing the following recursive relationships at (time) 
step n: 




(75) 


n " La* a* J (As) n 


(76) 




X_ + Y. 


(77) 


where: 


x„ = x. 


e -0.0455(.-M 2 )(A«V +0.165 (a n -a n _,) 


(78a) 


Y n ) e -o.3(i-M z )(A*)n +o.335(a n -a n _,) 


(78b) 


2Un 

(As)n 5 He 


(79) 


Here A^is azimuthal stepsize, Q is rotor speed, c is chord length, and 
is streamwise velocity. 
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The instantaneous angle-of -attack, a , is described in the section 

normal to the midchord, 6 and <f> being the pitch angle and inflow angle, 

n n 

respectively. The numerical calculation of the nondiraensional angle-of- 
attack rate. A, poses special problems. The nondimens ional time derivative of 
pitch angle in Equation (76), 90 n / 3^, may be computed analytically from 

the known control angle and elastic torsion response rates, whereas the time 
derivative of <f>must be computed using some form of numerical differentiation. 
The nominal method suggested in Reference 12 is a backward difference scheme. 
However, in some applications, this method was found to give violent numerical 
instabilities and an alternate method was required. The alternate method 
selected is based on the assumption of a predominantly oscillatory response 
at some user selected frequency, CJ, which typically would be taken as the 
dominant blade torsion natural frequency. These two numerical results are 
given below: 




{ 


— (I.5*n- 2* n _, + .5<£ n _ 2 ). w = 0 


<u cos oi AV' 
sin uJ AV> 


(cos w A'/' • <#>„- <#> n _ | ).< J j > 0 


(80) 


Prediction of Dynamic Stall Events 

In the present method it is considered important to accurately 
predict three major events associated with dynamic stall. These events, 
as shown in Figure 9b, are the stall onset, the vortex at the trailing edge, 
and the reattachment. The next section describes the semi-empirical equations 
that are used to predict these events. 

Onset of Stall 


Because the dynamic stall airloads acting on an airfoil are highly 
influenced by the leading edge vortex, an accurate prediction of the instant 
the vortex breaks away from the leading edge (moment stall point) becomes 
very important. The occurrence of moment stall depends on factors such as 
Mach number, the airfoil shape and the pitch rate. 

Under the conventional quasi-static theory formulation, the stall 
is assumed to occur when the effective angle-of-attack reaches the static 
stall angle. 


Em 


tt 


(81) 
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In general, a is assumed to vary with the airfoil shape, Mach number and 
Reynolds numblr. To some extent, the value of also depends on the 
criterion followed for stall. 

Under the present formulation, the relationship represented by 

Equation (81) is extended to include dynamic stall effects, and an assumption 

is made that at the dynamic stall point, in general, the effective angle of 

attack, Qt is not only a function of a , but also depends on the pitch 

Em ss 

rate at stall, A , and the instantaneous angle-of-attack at stall. That is, 
m 


a Em = F ^ a ss ’ A m’ a Drr? 


(82) 


The actual functionality F depends on the type of stall and on the criterion 
followed for stall. It is assumed that F varies with airfoil shape, Mach 
number, and Reynolds number, and can be established empirically. Linearization 
of the relationship of Equation (82) with respect to parameters and 

around quasi-steady conditions, o (1+e), leads to the following simple 

s s 

expression for o t , the angle at which dynamic moment stall first occurs: 

Dm 


W (lte +c »m A m + c .™ 0 .J 0 M 


(83) 


Here, a represents the value of the parameter, a , at the point of 
w w 

m 

moment stall. Thus, instead of the function F, one can determine empirically 

the coefficients e, C , and C for various Mach numbers, Reynolds numbers. 

Am wm 

and airfoils. In Equation (83), the last two terms represent the delay in 
dynamic stall when compared with quasi-static stall. Other available methods 
(References 15, 16) represent this delay in stall by a constant time delay. 
However, Equation (83) is a much more general relationship which predicts the 
onset of dynamic stall quite accurately for airfoils experiencing unsteady 
motion. 

Vortex at Trailing Edge 

Normally, after the occurrence of moment stall, there is a significant 
increase in negative pitching moment due to the travel of the stall vortex. 
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The maximum negative pitching moment occurs when the vortex is near the 
trailing edge of the airfoil. For the case shown in Figure 9b the 
instant when the vortex leaves the trailing edge is marked by ’ TE*. 
Preliminary results have led to the following empirical relationship for 
predicting the instant when the vortex leaves the airfoil: 


Smt = hO/(C A f Apm + Gpm) 


( 84 ) 


Here s is the total nondimens ional time for the vortex to travel from 
mt 

the leading edge to the trailing edge. Once again, the coefficients 
C and C vary with Mach number, airfoil shape, sweep, and Reynolds number, 
“t 

Reattachment 


The instant when the reattachment of the flow occurs is marked in 

Figure 9b. Normally, for low Mach numbers (M^0.4) the reattachment occurs 

at an angle a which is less than the static stall angle. At higher Mach 
RE 


numbers, where the static stall may be induced by shocks, the reattachment 

angle a can be higher than the static stall angle, a . In the present 
RE ss 

formulation, a general expression for a is assumed and is given by: 

RE 


a RE s 0 € + C AR A[) m + c WR a WM ) a ss 


(85) 


In general, for a given airfoil, the values of and C^, as used in 

Equation (85) for reattachment, are quite different from C and C used 

Am wm 

for stall onset. However, the value of the parameter e is the same in both 
of these equations. 

This completes the description of all the events associated with dynamic 
stall that are required to compute the unsteady stall aerodynamic characteris- 
tics of an airfoil. It should be noted that the present formulation does not 
require explicit prediction of so-called ’dynamic lift stall 1 . Normally a 
sudden loss of lift occurs due to increase in the relative distance between 
the stall vortex and the airfoil surface. These effects are included 
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implicitly in the formulation of the unsteady lift coefficient which is 
described next. 


Unsteady Section Coefficients 
Unsteady Lift Coefficient 

The unsteady lift coefficient. C . of an airfoil in the time domain 

Lu 

under the present synthes izat ion is described by the following expressions: 


Ci u = Crjta - Aa, - Ac^) + o OL Aci| + AC L i+ AC L 2 


( 86 ) 


Aa, = (P, A + P 2 a w + P 3 ) a 


ss 


( 87 ) 


Aa 2 = S 2 a ss 


( 88 ) 


AC L | = 0| A + o 2 a w ♦ o 3 (a /a ss ) + Q 4 (a / a^) 


( 89 ) 


ac l2 =0 5 8, +• o 6 Aa 2 + q 7 (a Dm ) 




(^i s m) 


( 90 ) 
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(91) 


s m * 


zult-W 


c 


0 

a < ass 

(a/a #s - 1) 

a ss < a < a om 

(<W a *s -«> 

[l.O-(s m /s mt ) J 

0 

s m > s mt 


(92) 


(a/a ss -D 


a -a 


^ a on/ a *s - 1 ^ a TE -a 


RE 


a - a ss 

a « Z a * a om 

0<S m <S m1 


a RE - Q - a TE 


a < a 


RE 


(93) 


The synthesized unsteady lift coefficient (Equation (86)) has been 
expressed as a sum of static at some shifted angle (a-Aa^ ~ 

plus an incremental lift coefficient ( AC^ + AC^) • an Sl e 

is given by Equations (87) and (88) and the incremental lift coefficient 
by Equations (89) and (90). (The quantity a in Equation (86) is the 

OL 

conventional static lift curve slope.) The shift in angle (Equation (87)) 

is present even when no stall occurs, and the A a shift in angle 

(Equation (88)) is mainly associated with the occurrence of dynamic stall and 
subsequent reattachment. Similarly, the A C (Equation (89)) represents 

essentially the unsteady effects over static C for dynamically unstalled 

LD 

airfoils, and A C (Equation (90)) represents the effects associated with 
LZ 

the dynamic stall events such as vortex formation and reattachment. In fact, 
the last term in Equation (90) represents explicitly the suction effects of 
the leading edge vortex and equals zero when no vortex exists. 
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Thus, Equation (86) is a general expression for unsteady C 


L 


even when no 


dynamic stall occurs. For unstalled cases, 
A 0^2 are essentially zero. 


the magnitudes of A a 


2 


and 


The 6 ^ parameter in Equation (90) is an empirically determined 

constant and is nominally equal to 0.18, The quantity s , as described 

m 

by Equation (91), represents the nondimens ional time measured from the 
instant of the occurrence of dynamic moment stall. The unknown parameters 
through and through are determined empirically by means of a 

least-squares curve-fitting of Equation (86) with the test data. It should 

be noted that most of the terms in Equation (86) are linear in parameters ot» 

A, and a . 

w 


Unsteady Moment Coefficient 


The unsteady pitching moment coefficient, C , has been formulated 

Mu 

using relationships similar to those for C and is described below: 

Lu 


C Mu = C M$^ a “ Aa z ) ♦ °om + AC m 


(94) 


AC m s ij, a + i , 2 a w + i, 3 (a/a $# ) + i , 4 |a w | 

+ i, 5 8, +i, 6 ^+1,^^^ 


Here a represents the static pitching moment slope at zero angle-of- 
om 

attack and it normally equals zero. The last term in Equation (95) represents 
the vortex effects. For unstalled airfoils, the last three terms in 
Equation (95) are zero. The unknown parameters through once more are 

determined by the least-square curve-fitting of Equation (94) to the 
test data. 
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Unsteady Drag Coefficient 


The unsteady drag coefficient, C , appears to vary with the dynamic 

Du 

parameters in the same way as C and is described as follows: 

Mu 


^ou * Cqs ~ AGtg) + AC q 


(96) 


AC 0 =R, a + r 2 a w + R 3 (a/a ss ) + R 4 |a w | 

+ R 5 8 , + R 6 S 4 + R 7 Aa 2 + Rg a Dm A 0 m s m 


where : 


83 


!° a * a ss 

(a/a„-i) OsgSasao,,, 

( a om^ a ss”*^ ® - s m - s mt 

0 S m >S mt 


(97) 


(98) 


8 «=< 


(a/a ss - 1 ) 


toon/ass-') 


a — ®ss 
a„< a< a D m 

[*• — ^ ^ ® — s m — ^ mt 


0 


> ^ mt 


(99) 


The last term in Equation (97) represents the effects of the stall 

vortex on the unsteady drag. For unstalled conditions, the last four 

terms in Equation (97) are essentially equal to zero. Once more, the 

unknown parameters R through R are computed using the linear least-squares 

1 8 

curve-fitting of Equation (97) to the unsteady drag test data. 
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Because the higher Mach number unsteady data used to synthesize the 
coefficients in the G400PROP code did not include unsteady drag, the 
use of unsteady drag is not available to this version of G400. 

The description of unsteady drag coefficient given above was included 
herein for completeness. 

Description of Additional Synthesization 

The empirical relationships, for the prediction of stall events 
(Equations (83) through (85)) and for the description of unsteady airfoil 
characteristics (Equations (86) through (99)), have been established by 
utilizing a large number of available oscillatory airfoil test data sets. 
Furthermore, by illustrating the excellent correlation between the test 
and synthesized results, the generality of these empirical relationships 
to adequately represent the effects of variations in Mach number, sweep, 
and airfoil shape has been clearly demonstrated (References 12 and 17). 

This section describes the similar correlation results obtained under the 
present study, which relates mainly to the synthesization of the high 
subsonic Mach number data. 

Test Data Used for Present Synthesis 

The first step in the procedure for synthesis normally involves 
preparing a data set consisting of the loop data obtained for an airfoil 
at the same Mach number, Reynolds number, and sweep angle. Normally, a 
set of fifteen loops, consisting of both unstalled and stalled data, is 
found to be sufficient to establish the values of the empirical coefficients. 
The second step of the synthesis procedure consists of determining the 
empirical coefficients through least-squares fitting. The final step 
involves reconstructing the data from the empirical relations and comparing 
the synthesized data with test data. 

Table I provides a list of all the data sets that were successfully 
synthesized under the present study. The data sets listed in Table I were 
acquired from two different sources (1) NASA CR-2915 (Reference 13) containing 
data sets 1 through 3, and (2) USAAVLABS TR-68-13B (Reference 14) containing 
data sets 4 through 7. 

Each of the seven data sets represents a unique combination of test 
conditions. As a result, the values of the various empirical coefficients 
obtained are, in general, different for each of these data sets. Also, it 
should be noted that each of these data sets have, in general, a different 
static airfoil characteristic associated with them (steady state C^, C^, 

Cp variation with a) . 
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TABLE I: TEST DATA SETS USED IN AIRFOIL UNSTEADY STALL SYNTHESIZATION 


Data 

Set 

No, 

Source 
Ref. No. 

Airfoil 

Type 

Mach 

No. 

Reynolds No. 
xlO -6 

Parameter Range of Test 

Data Used 

k 

a 

o 

a 

1 

13 

NLR-1 

0.5 

6.3 

0.07-0.22 

0.0-12.5 

2. 5-7. 5 

2 

13 

NLR-1 

0.6 

9.1 

0.06-0.18 

0.0-12.5 

2. 5-7. 5 

3 

13 

NLR-1 

0.7 

10.1 

0.05-0.16 

0.0-12.5 

2. 5-7. 5 

4 

14 

V0012 

0.4 

4.8 

0.0-0.31 

5.0-15.0 

2.5-7. 5 

5 

14 

V0012 

0.6 

6.2 

0.0-0.25 

5.0-10.0 

2. 5-7. 5 

6 

14 

V2301-158 

0.4 

4.8 

0.0-0.25 

5.0-15.0 

2. 5-7. 5 

7 

14 

V2301-158 

0.6 

6.2 

0.0-0.25 

5.0-10.0 

2. 5-7. 5 
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Comparison of Synthesized Loop Data With Test Data 


This subsection discussed the results obtained from the curve-fitting 
of Equations ( 86 ) and (94) to the test loop data corresponding to lift 
coefficient and pitching moment coefficient, respectively. As a typical 
case, consider all the lift coefficient loop data contained in data set 
number 3 (see Table I). When these loop data are curve fitted to 
Equation ( 86 ) in a least-squares sense, the values of unknown parameters 


through 
P^ through P^ 


and Q through 
and Q 1 through 


are obtained. When the values of the 
parameters are inserted in Equations ( 86 ) 


the resulting time domain equation represents the two-dimensional unsteady 
lift coefficient of the NLR-1 airfoil at Mach number 0.7 for essentially 
all dynamic conditions. 


To illustrate the accuracy of the resulting equation, a sample of the 
loop data for this case has been reconstructed from the equations and the 
comparisons of these synthesized loops with test data are shown in 

Figure 10a, The differences between the test data and the synthesized are 
small and these differences are comparable to test data accuracy. 


Similarly, when all the pitching moment coefficient loops contained 
in the data set number 3 in Table I are curve-fitted to Equation (94), 
the values of unknown parameters through n 7 are obtained. The comparison 
of the synthesized C loops with test data is shown in Figure 10b. Once 

" Mu 

again, the reconstructed loops match very well with the test data. The 

maximum negative C is generally predicted accurately for all the stalled 
Mu 

loops. 

Similar computations for the six other data sets contained in Table I 
have been successfully carried out. Figures 10a through 12b illustrate the 
good agreement obtained between the synthesized loop data and the test data 
corresponding to the highest Mach number data sets for each of the three 
airfoils. These figures correspond to the published data contained in 
References 13 and 14 as obtained in the Boeing two-dimensional, variable- 
density wind tunnel. It should be noted that Figures 10 through 12 present 
only lift and pitching moment loops because no drag data were included 
in these references. 
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DATA, COMPUTED 



Figure 10a. Comparison of Synthesized Lift Coefficient Loops with Test Data, NLR-1 Airfoil, 
M =0.7, R n = 10.0 x 10 6 Data Set No. 3 in Table I 





COMPUTED 



Figure 10b. Comparison of Synthesized Pitching Moment Coefficient Loops with Test Data, 
NLR-1 Airfoil, M =0.7, R n = 1.0 x 10 6 Data Set No. 3 in Table I 








Computed 
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Figure 11a. Comparison of Sythesized Lift Coefficient Loops with Test Data Vertol Modified 
NACA 0012 Airfoil, M =0.60, R n =6.2 x 10 6 Data Set No. 5 in Table I 






Data, Computed 



Figure 11b. Comparison of Synthesized Pitching Moment Coefficient Loops with Test Data Vertoi 
Modified NACA 0012 Airfoil, M =0.60, R n = 6.2 x 10 6 Data Set No. 5 in Table I 







Data, Computed 



73 


ALPHA, DEG. 

Figure 12a. Comparison of Synthesized Lift Coefficient Loops with Test Data V2301-1.58 Airfoil, 
M =0.60, R n =6.2 x 10 6 Data Set No. 7 in Table I 




Data, Computed 


i 
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Figure 12b. Comparison of Sythesized Pitching Moment Coefficient Loops with Test Data 
V2301-1.58 Airfoil, M =0.60, R n = 6.2 x 10 6 Data Set No. 7 in Table I 


UNSTEADY UNSTALLED SUBSONIC AERODYNAMICS 


The use of quasi-static airloads in the aeroelastic analysis of 
advanced propeller systems (such as prop-fans) lends itself to simplicity and, 
hence, economy rather than accuracy. For an accurate quantitative aero- 
elastic analysis, unsteady aerodynamic forces become indispensible. This 
can be seen by noting the lift coefficient variations with reduced frequency 
shown in Fig. 13 for a two-dimensional airfoil at a subsonic Mach number 
typical of prop-fan operations. The reduced frequency range shown in Fig. 13, 
moreover, is typical of the vibration modes of real prop-fan blades. The 
aerodynamic force lag is substantial as implied by the imaginary part of the 
lift coeffiqient. 

The majority of the available unsteady aerodynamic lift and moment in- 
formation for airfoils comes from theory or experiment in the (real) fre- 
quency domain instead of in the time domain. This is mostly due to the sim- 
plicity in mathematics and experimental effort in working in the frequency 
domain. In order to perform time-history solutions for an aeroelastic prob- 
lem, however, the frequency domain unsteady aerodynamic data must be properly 
transformed into the time domain. The frequency domain unsteady aerodynamic 
data are typically in tabulated or transcendental function form. As a 
result, it is difficult to perform a transformation which is both accurate 
and economical. In cases where only eigensolutions are required, there 
remains a fundamental problem of generalizing data available only in the fre- 
quency domain (constant amplitude oscillation) to the complex frequency or 
Laplace variable domain (decaying and growing oscillations) . 

In order to overcome the above mentioned difficulties, Pade approximants 
have been introduced in the literature as an approximate but consistent way 
to bridge the gap between the (real) frequency domain usnteady aerodynamic 
data and the time domain description of the unsteady aerodynamic forces. See, 
for example. Reference 19. As opposed to the generally transcendental nature 
of the unsteady aerodynamic data, the Pade approximants are defined in terms 
of rational functions that are known to have simple Laplace inversions or 
inverse Fourier transforms. Besides its mathematical advantage, the Pade 
approximant also provides a quick method for interpolating and/or extrapolat- 
ing the frequency domain data, which are usually limited to some discrete 
frequencies . 

The sources of unsteady aerodynamic data for generating Pade approxi- 
mants can be theoretical and/or experimental in nature. If necessary, the 
data source can even be nonlinear as exemplified by the time domain transonic 
LTRAN2 code or its advanced versions (References 20 and 21). Such nonlinear 
sources can be used either in the frequency domain by explicitly making the 
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UNSTEADY LIFT COEFFICIENTS C fl 





Figure 13. Example of Unsteady Aerodynamic Airfoil Characteristics and 
Associated Pad& Approximation 
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airfoil motion sinusoidal (in which case only the first harmonic component 
of the total airfoil is extracted and used) , or in the time domain itself 
using indicial responses* In the latter case, the Pade form lends itself 
equally well to an exponential fitting procedure such as Prony's method 
given in Reference 22. 

In the following subsections, the sources of unsteady aerodynamic data 
used in this study are first described and then the data synthesization pro- 
cedures for rendering these data to Pade forms are discussed. Then, in the 
subsequent subsections, the details of going from the Pade forms to linear 
differential equations are described. 


Sources of Unsteady Airloads 

» 

Either theoretical or experimental data can be used in Pade approxima- 
tions. In this report, both data source types were used. The theoretical 
linear unsteady aerodynamic data source selected is the work of Jordan 
(Reference 18) for two dimensional flow about an isolated flat plate air- 
foil. The experimental data source used is that of Davis and Malcolm 
(Reference 23) for the NACA 64A10 airfoil. The data were put in standard 
forms according to the following lift and (quarter chord) pitching moment 
coefficient definitions before synthesization. 


L 

^/>u 2 h 


(100a) 


l-/ou 2 ca 


M 


'nn k - 


^/>u 2 ch 


(100b) 


(100c) 


Cma = i./>u 2 c 2 a (100d) 

Here h and a are the airfoil plunging amplitude and the airfoil pitching 
amplitude, respectively. The pitching motion and the moment are defined 
throughout this section about the quarter chord. The total lift and moment 
coefficients are then given simply by the sums of the plunging and pitching 
results. 

* 
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c i = C Xh h + c ia a 


(101) 


Cm = C mh h + C ma a 


( 102 ) 


Synthesization of Data to Pade Form 

! 

The Pade approximant of any unsteady frequency domain aerodynamic 
coefficient C is defined as the ratio of two polynomials in complex fre- 
q uenc y . N+( 

j=i 

c(w) s — ^ : do3) 

y b, (i^) J 

A J 

where N is the order of the approximant. The degree of the denominator is 
lower than that of the numerator by one, because of the known asymptotic 
behavior at large frequencies. 

Physical Constraints 

The Pade coefficients aj and bj in Eq . (103) are determined by imposing 
the following requirements: 

(a) The zero frequency data must be satisfied exactly. 

(b) The Pade approximants should approach the piston theory results 
asymptotically for large frequencies. 

(c) The available data (except for zero frequency) will be approximated 
by Eq. (103) in the least-square sense. 

(d) The resultant poles must be stable. 

Using basic concepts given in Reference 24, the lift and moment coeffi- 
cients of plunging and pitching airfoils based on the piston theory can be 
written in the following form: 


C = A + 



B 


(104) 
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where A and B are frequency independent. These constants (A, B) for both 
lift and moment coefficients of plunging and pitching airfoils are shown 
below. 


Constants in Force and Moment Coefficients 
From Piston Theory 


Unsteady 

Coefficient 

A 

B 

c *h 

0 

4 

M 

Cm h 

0 

1(1 - h) 

M ' 2 C ' 

Zla 

4 

M 

4.1 x 0 
M ' 2 C > 

r 

4 , x s 1 v 

4,1 x 0 1 x s 1 x s X 0 


M ~ 

M ( Z~~ + 2 T’ 3 " T 


Note: x Q = pitching axis location measured from L.E. 

x g = point about which moment is taken 
c = full chord 


Synthesization Techniques - Least-Square Fit and Weighting Schemes 

i 

For simplicity, a two pole Pade approximant is sought, namely, N = 2. 


_ . o(iou) 3 + b(ia>) 2 + c(M+d 

C(cu)= 

(iu>) 2 + e(io>) + f 


(105) 


To satisfy the zero frequency and high frequency limits we have, for x g = 
x Q = c/4, the following constraints: 

d = c ( 0 ) f 

0= i//3 


where : 

( -M/4 for Ci h 

£= | M tor C ia 

I M for Cm h 

' — I2M/7 for Cm a 


(106) 


79 



Therefore, we are left with four unknown constants: b, c, e and f. Let 


C (oi ) = Cr (cu) + i Ci(w) 


and multiply out Eq. (105) to yield the following real equations. 


cu 2 0 -cuCj C r -C( 0) 
0 cu cuC R - C 2 



(107) 


I± e il s JL 

The unsteady aerodynamic coefficient C(co) is assumed to be known for n 
frequency values. Equation (107) in general cannot be satisfied for all fre- 
quencies when n > 2. In fact, we would have the following 2n equations 
for for unknowns: 


[A] M | B f 

2nx4 4xi znxi 


(108a) 


where: 


w 



(108b) 


The least-square solution for {x}, however, can be found as the solution of 
the set of modified equations. 


[E] |«| = |0| 

4x4 4x1 4x1 


(109) 


where: 

[E]= [A] T [A] ; [D] = [A] T [D] 


» 

Equation (109) is the final required form used to generate the Pade 
coefficients for both the Jordan theoretical data and the Davis and Malcolm 
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experimental data, as supplemented by the Jordan data for missing 
frequencies. It was initially found that the Pade approximants of some 
of the lift coefficients would not provide stable poles and almost all 
Pade approximants of the moment coefficients contained unstable poles. 

^ightinj^ Schemes_ 

It was believed that the unstable poles might have been caused by the 
undue weighting to high frequency data in the least-square procedure de- 
fined by Eq. (109). This can be seen in the coefficients in Eq. (107). 
Therefore, several attempts to minimize this undue weighting for high fre- 
quencies were made by dividing the two real equations in Eq. (107) by 
several chosen functions of as follows: 


Option 

Divide 1st Equation by 

Divide 2nd Equation by 


2 


1 

0) 

a; 

2 

2 t 

u) In a) 

w In a) 

3 

(a) In a) 

a) In a) 

4 

only for ca > 1 

a) only for to > 1 

5 

o)^ only for w > 1 

only for oa > 1 


As a result of these normalization processes, almost all lift 
coefficients and most moment coefficients resulted in stable Pade approxi- 
mants. A further interpolation and extrapolation procedure applied to the 
poles rendered all Pade approximants stable as required. 


? 

Working Forms of Pade Approximants - Partial Fractions 


The Pade approximant in Eq. (103) can be written in terms of its par- 
tial fractions. Then the total lift coefficient becomes, according to 


Eq. (101), # 

C£= (A, + 


A 2 P 

P + P, 


A 3 p b 3 p 

+ ~~T )++ (B.p+B 2 + -f- 

P"P 2 P-P, 



( 110 ) 


where we have formally replaced the Fourier transform variable iai by the 
Laplace transform variable p. Meanwhile, the inflow angle, <f>, replaced the 
plunging velocity variable, ph/c, and the pitch angle, 0, has replaced old 
symbol a. A similar equation holds for the total moment coefficient with 
generally different constant coefficients and poles in Eq. (9). These 
coefficients and poles are summarized in Tables IIA through IID. 
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TABLE IIA - PADE APPROXIMANT COEFFICIENTS FOR C 


JC 
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TABLE IIB - PADE APPROXIMANT COEFFICIENTS FOR C 
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TABLE IIC - PADE APPROXIMANT COEFFICIENTS FOR Cj, 
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TABLE I ID - PADE APPROXIMANT COEFFICIENTS FOR C 
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Relation of Pade Approximation to Unsteady Decay Parameter 

The introduction of the formalized Pade form to the description of 
unsteady airloads is but a generalization of an approximation to the Wagner 
function originally formulated by R. T. Jones (Reference 25). In the pre- 
ceeding section, the unsteady decay parameter, a w , was defined using a 
generalized form of this approximation (see Eq. (68)). The resulting 
’’effective 11 angle-of-attack, otg, defined using this unsteady decay parameter, 
is a useful aerodynamic tool in its own right. It can be used independent 
of the unsteady stalled airloads theory to approximate low frequency un- 
stalled unsteady airloads. Although the effective angle-of-attack concept 
assumes that plunging motion can be treated as an equivalent pitch angle, 

Eq. (66) can still be used to formulate an airloads description similar to 
that given above in Eqs. (101) and (102): 

c i= f (•§•)« + c ia„°E <m> 


Cm c/« ' ' T 1 ® + (Cm ao ) c/4 % 


( 112 ) 


It can be shown that the Laplace transforms of Eqs. (Ill) and (112) 
together with that for the effective angle-of-attack, oig, result in mathe- 
matical forms which are identical to Eq. (110). The various constants de- 
fining the Laplace transformed lift and moment equations, for the Generalized 
Wagner function (oig) formulation are given in Table III. 


Differential Equation Form 

The starting points for formulating practical differential equations 
for the airloads are the Laplace transformed equations for lift and pitching 
moment, as typified by Eq . (110). The expressions for both the lift and 
moment coefficients may then be rewritten in the same abbreviated general 
form: 


AC = [a,<#> + A 2 X+A 3 y + B,p0 + B z 6 + B 3 


(113) 
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Table III. Pad4 Coefficients for Generalized Wagner Function 


Symbol Cj^ C 


-0.022758 2 

-0.022758 2 

I 

-0.156 2 

-0.158 2 

C ^a 

o 

C i"ao 

-Co (0.165) 

o 

-C ma (0.165) 

-Cp (0.335) 

x a 

o 

-C^ (0.335) 
o 

C £ 

Cm 

a ! 

-0.022756 2 

-0.022758 2 

-0.156 2 

-0.158 2 

Cp +0.5:: 
a 

o 

C - 0.25ti 

ma 

o 

C £ 

a 

° 2 
-Cp (. 1648-. 00758 ) 

‘■a 

° 2 
-C { (.33517 - .100558 ) 

o 

Cm a 

o 

-C_ (.16483 - .00756 2 ) 

“a 

0 2 
-C_ (. 33517-. 100558 ) 

“a 

o 





where Ac refers to a perturbational implementation of the Pade theory, to 
be discussed in a later section, and where the augmented state variables 
appearing in Eq . (113) are given by: 


x = 

9<t> 

P-* 

(114) 

y = 

p<*> 

A 

(115) 

z = 

P-P 2 

99 

(116) 

w = 

p-R 

P0 

(117) 


P-P 2 


The above constants and poles have unique values for lift and for moment 
as shown in Tables II and III. 

For each coefficient, there are generally four Laplace transformed 
expressions to solve of the form presented in the above equations: 


P4> 



Since the Laplace operator is invertible, the associated differential 
equation becomes: 


Px- P<£ = 


A 

P ( x 


(118) 


and then after nondimensionalizing the time differentials by aerodynamic 
time (chord/velocity) and rearranging: 



(119a) 


where U and c are the nondimensionalized velocity and chord, respectively. 
Similarly: 
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(119b) 


HIrW 


(119c) 

(119d) 


Solving these differential equations over a time step, Aip, gives the 
following formulas: 


An.. 

R ■= A* 


x k = x k _, e 


up,/c 


* A U , 


(120a) 
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W k= w k-. e C 




+ — 


* 

e 


up 2 /c 


[.***♦-.] 


(120d) 


where (k-1) and (k) refer to successive time steps, where AiJ) is the time 
increment, and where each expression is formulated for both lift and for 
moment. 


Implementation Within the Time-History Solution 

i 

The Pade coefficients and poles, computed using the data synthesiza- 
tion techniques discussed above, were included in the G400PROP time-history 
solution. The Pade approximations of the Jordan (theoretical) and the 
Davis and Malcolm (experimental) aerodynamic data sets were tabulated for 
variations in Mach number from 0.5 to 0.95 for Jordan, and for a Mach number 
of 0.8 for Davis and Malcolm. A linear interpolation scheme was devised in 
order to obtain Pade values at intermediary Mach numbers. 
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Since the Pade theory is based upon classical small amplitude 
aerodynamic theory, a perturbational implementation appears warranted. To 
this end, an option was included to use either perturbational or total 
(perturbational plus static) pitch and inflow angles in the calculations 
of Ac£ and Ac m . The option to use the total angle description in place of 
the perturbational angle description was included to allow for lack of 
static airfoil data. 

With the perturbational approach, the steady-state values of the pitch 
and inflow angles are estimated by eliminating the contributions to these 
angles due to elastic responses. This estimation is accomplished by elim- 
inating the elastic torsion from the pitch angle and the explicit rate 
dependent terms from the tangential and perpendicular section velocities 
used to form the inflow angle. Thus, referring to Eqs . (54), (46) and 
(50) we have: 

NTM 

j=l 

6 

Therefore the perturbational angles become: 

e = e- 9 0 ; * s * - 


(121a) 


(121b) 


at each time step. 

These essentially filtered values of 6 and $ would then be used in 
computing the Pade Ac^ and Ac coefficients, and the 0 and $ values would 
be used to compute a steady-state c^ and c m (from the static airfoil 
tables). The static c^ and c^ wou?d be added to the Pade Ac^ and Ac m : 

Cj £ = + AC^ } Cm ” Cm 0 + ACm (123) 


The object of this filtering is to produce perturbational angles for 
use in the Pade calculations which calculate only the transient airloads. 
Therefore, the filtered option should produce lift and moment coefficients 
which, with time, approach those of the quasi-steady case. The above 
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approach therefore uses only perturbatlonal angles-of-attack in computing 
the unsteady airloads together with estimated steady-state angles-of-attack 
in computing the steady airloads to determine the overall transient aero- 
elastic responses. 
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COUPLING WITH PROPELLER/NACELLE PERFORMANCE ANALYSIS (PANPER) 


Overview and Interfacing Issues 

PANPER is an analysis developed at UTRC for high speed propeller- 
nacelle aerodynamic performance prediction*. The program produces in- 
duced velocities, in cylindrical coordinates, for use in blade response 
programs such as G400PROP . PANPER assumes axisymmetric flow and single 
propeller operation, as opposed to coaxial, counter-rotation duel propeller 
operation. The induced velocity field provided by PANPER includes the 
flow distortion effects of the nacelle. 

In usage, the G400PROP analysis must be run first, to obtain the 
various geometric and aerodynamic parameters required by the variable in- 
flow analysis. Since PANPER assumes axisymmetry, these parameters must be 
the azimuthally averaged values evaluated at the last rotor revolution. 

A transformation must be made between the ”5” coordinate system (preconed 
and prelead-lagged feathering axis) and the rectilinear coordinate system 
required by PANPER. Also, the airfoil data used in G400PR0P must be 
modeled in polynomial format for use in PANPER. 


Geometric and Aerodynamic Data Required for PANPER 

Nondistributed Geometric and Aerodynamic Descriptors 

The principal nondistributed parameters which must be conveyed from 
G400PROP to PANPER are as follows: 

1. The number of blades. 

2. The blade radius, R, feet. 

3. The speed of sound, a^, feet/second. 

4. The density of air, p, lb-sec^/ft^. 

5. The forward flight speed, Vj, knots. 

6. The propeller rotational speed, ft, rpm. 

7. The momentum induced velocity, v im , defined positive in the upflow 
direction, feet/second. 


This program was developed under Contract NAS3-20961 and is to be docu- 
mented in a two volume set entitled "An Analysis for High Speed Propeller- 
Nacelle Aerodynamic Performance Prediction". 
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Distributed Geometric Quantities 


The PANPER routine requires distributions of the spanwise segment 
centers and boundary points at the quarter chord, and in its rectilinear 
coordinate system. It also requires distributions of the chord and pitch 
angle in the same coordinate system, referred to herein as the "4" system. 


The following material draws upon the earlier sections relating to 
structural and aerodynamic sweep, and refers to the details shown in Fig. 

14. As shown in this figure, the "4" coordinate system is defined by vec- 
tors taken locally at the intersection of the segment midchord with the seg- 
ment boundary. L is the local radius of the midchord-boundary point P, 
measured from the hub; £ is the angle which this radius makes with the pre- 
coned and prelead-lagged feathering axis. The distance measured from the 
swept G400PROP segment boundaries to the "4 M system segment boundaries is 
referred to as n* where n varies along the chord from 0 at the midchord to 
n LE at l ea< H n 8 edge and at the trailing edge. 


First, the total "5" system deflections must be formed, based upon a 
combination of the elastic axis location and the elastic bending terms as 
described in Eqs . (27) and (28), and upon a transformation due to aero- 
dynamic sweep of the arbitrary chordwise locations within the airfoil sec- 
tion. As presented in an above section (see Eq. (44)), the aerodynamic 
sweep transformation matrix can be defined as: 



where: 
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(A) 


(A) 
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-smAfsX 
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• * (A) * V 

smAf sinAgi 


(A) 


cos A 


(A) 


, . (A) 

smAfs 


. (A) 

cos Af 5 


(124) 


,<A) 




sin 2 A^— sin 2 A' 


(A) 

h 


(125) 


and where the aerodynamic sweep angles are as defined by Eqs. (42) and (43). 
The arbitrary chordwise locations within the airfoil section are formed 
using straightforward trigonometric rotations. Using the leading edge of the 
airfoil section provides the following example: 
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Figure 14. Schematics of the “4” and “5” Coordinate Systems Relating to 
G400PROP/PANPER Interfacing 










(126) 


( X *) (*») ( 0 

r 5 t = V5 ( + [ TAS<A> 1 j y 'o L E cos0 

^ z 5 Ve V z 5 ) ( Yio le sin® 

where: 

y '0 L E =y i0 0C + C/4 < 127 > 

and where y^Oor* t ^ ie ec *gewise location of the section quarter chord point 
relative to tne elastic axis, is input to the program. Similarly, the 
total "5" system displacements may be found at other chordwise locations 
within the airfoil section such as the midchord, the quarter chord, the 
trailing edge, etc. 

Next, in the "5" system, the intersection of the section midchord with 
the section boundary must be obtained for use in determining the angle £ 
and the distance L. The inboard boundary displacements are assumed zero. 
Direct interpolation is required to obtain the midchord boundary positions 
at arbitrary spanwise locations (n) along the span, and extrapolation 
must be used to obtain the outermost midchord boundary point. 


Zs MCB^ " 


Z5 Mc^ n ^ Ax(n- 1) 
Ax(n) + Ax(n- 1 ) 


z 5 MC ( n -^ 'Ax(n) 
Ax(n) + Ax(n- l) 


(128) 


Zs MCB^ Z 5 mC^ + ( Z 5MC^ Z 5MCB^ N 


(129) 


where MCB refers to the intersection of the segment midchord with the seg- 
ment boundary, MC refers to the intersection of the segment midchord with 
the segment center, the parenthetical indices refer to the segment number, 
and (N) refers to the last spanwise segment. 

Solving for the angle £ and the distance L: 

, . r W n) - esin8 » 

t = ton 

L ecos8 B + x 5 (n)cos/3 B - z s (n )/3 B 



L = V(y 5MCB (n) - esin8 B ) 2 + (ecos8 B + X 5 MCB ^ cos£ B - Zsmcb^^b^ 2 


(131) 


95 



where Fig. 14 displays the various geometric quantities such as offset* e, 
prelead-lag, 6 fi , precone, B fi , etc. 

Once £ and L have been defined, the "5" system segment boundary dis- 
placements may be found, by interpolation, at points along the chord other 
than the midchord (n MC ** 0). For example, the segment boundary points at 
the leading edge of the chord may be found based upon the segment midpoints 
at the leading edge of the chord as follows: 



^LE + 


Ax(n) 


) 


LEB 


4-(Ax(n+i)+ Ax(n)) 


(n) 



(^LE - 


Ax(n+i) 


(n+i) 


Y (Ax(n+l)+ Ax(n)) 


(132) 


‘le 


(n) 


where LEB refers to the intersection of the segment leading edge with the 
segment boundary, LE refers to the intersection of the segment leading edge 
with the segment center, and n and n + 1 represent successive spanwise seg- 
ments. The center leading edge displacements were calculated in Eq. (126). 
However, remains so far unknown, the solution for which needs the 

development which follows. (A similar procedure might be followed for any 
other chordwise position within the airfoil section.) 

The coordinate transformation from the "5” system to the ”4" system is 
given by: 


(* 4 ) 

COS £ b COS £ 

sin i 

- /3 b cos£" 

/ x 5 + ecos8 B \ 
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- cos£ B sin£ 
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£ B sin£ 
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0 

cos/3 b 

( , ) 

z 5 


This transformation may be seen from Fig. 14 to be the product of a trans- 
formation from the "5" system to the hub coordinate system and a transforma- 
tion from the hub system to the rectilinear M 4 M system. 

The relationship = L, approximating the streamwise chord line 
through P, is given by the first row of Eq. (133). In general form: 


L = |cos£ b cos£, sin£, - 0 cos£]{x,+ *7X2} 


(134) 


where: 
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Equation (134) represents a linear relationship in n from which, for 
example, the values for and n T g may be determined. The interpolation 

procedure of Eq. (132) may then be completed. Having once determined the 
displacements of the segment boundary points in the "5" system, the boundary 
point values in the "4" system may be found using Eq. (133). 

The intersection of the segment quarter chord with the segment center, 
in the "4" coordinate system, may be determined to be midway between the 
intersections of the segment quarter chord with the segment boundary points. 
Finally, the chord and pitch angles in the "4" system are defined as follows 


2 -J Ka>- y -T E B )2 + (z 
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Distributed Aerodynamic Quantities 

The PANPER routine requires that lift and drag coefficients from 
G400PROP be in polynomial (quadratic) form about the operating angle-of- 
attack. Thus: 

C jT C Ao + C Ai a+C £2 aZ (138a) 
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(138b) 


C d = C d 0 +C d, a+C d 2 a2 

where a varies - 2° about the operating condition, and where the 
coefficients are determined using a Lagrange interpolation approach 
(Reference 22). It can be shown that: 
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(139a) 


(139b) 


(139c) 


where f(x^), f(x 2 ), and f(x^) are the lift coefficients computed at (a - 2°), 
a, and (a + 2°), respectively. Similarly, a series form of the drag coeffi- 
cient may be produced. 

The above lift and drag coefficients, along with the geometric quan- 
tities and distributions computed earlier, are written into a file chosen 
by the user for communication between G400PROP and PANPER. 


Assimilation of PANPER Generated Components of Variable Inflow 

The PANPER routine uses the geometric and aerodynamic data provided by 
G400PR0P to produce induced velocities, in cylindrical coordinates, due to 
the rotor, rotor wake, and nacelle. The radial induced velocity is defined 
positive from root to tip, the tangential velocity is defined positive 
opposite to the direction of rotation, and the axial velocity is defined 
positive in the normally thrusting direction. 

To convert the induced velocities back from the rectilinear "4" 
coordinate system to G400PROP f s "5" system, the inverse of the rotation 
matrix given in Eq. (133) may be applied. The new transformation matrix 
is defined as follows: 
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cos$ B cos^ me 

sin^mc 


-cos£ B sin£mc 

cos *mc 

S ' n 


^B 

0 

COS 0s 


(140) 


where £ MC is defined in a similar manner to £ , with the exception of its 
being taken at the midchord center rather than the midchord boundary. The 
above rotation varies with spanwise location, and is input to PANPER along 
with the geometric and aerodynamic distributions described above. PANPER 
uses this transformation to preconvert the calculated induced velocities 
to "5" system velocities at the G400PROP segment centers, before the data 
is transferred to G400PROP. These components of the variable inflow are 
then included in the airload calculations of both the eigensolution and 
the time-history solution. 
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EIGENSOLUTION 


Review of Basic Methodology 

All of the development presented in the previous sections has been 
directed to the accurate modeling of the various aeroelastic elements of 
rotary wings, and in particular, propellers. This modeling ultimately takes 
the form of a collection of nonlinear differential equations of motion. The 
two basic types of solution adopted by the G400 analyses for these differen- 
tial equations are eigensolutions and time-history solutions. 

Basic Solution Types 

The purposes of these two solution types are different, yet complementary. 
The purpose of the eigensolution is to calculate those inherent, general 
characteristics relating to vibration susceptibility (natural frequency and 
mode shape), and to stability (characteristic exponent and/or neutral 
stability point) . The eigensolution is essentially a solution for the so- 
called "homogeneous" solution to the differential equations. The purpose of 
the time-history solution, on the other hand, is to calculate the specific 
responses due to either self-excitation (together with appropriate initial 
conditions) , or environmental influences (control inputs and/or airflow 
distortions) . 

Of the two, the time-history solution is by far the easier in that all 
nonlinearities can be easily retained and a compact, implicit form of the 
equations can be used together with a variety of time-marching integration 
algorithms. The time-history solution selected for the G400 analyses 
(including G400PROP) is relatively simple and is well documented in 
Reference 1. Nothing substantially new has been added to this type of 
solution in the present study and further descriptions of time-history 
solutions per se are, therefore, omitted. 

General Characteristics of Eigensolution 


The principal reason eigensolutions are relatively more difficult to 
achieve is the requirement that the differential equations must be explicitly 
linearized. Hence, all the compactness afforded by the implicit formulations 
of nonlinearities available with time-history solutions must to a great extent 
be abandoned. The great advantage bought with the increased complexity of 
setting up an eigensolution, however, is the relatively low cost of solving 
the eigensolution compared with time-history solutions. 

Within the context of aeroelasticity , there are two basic types of 
eigensolutions: frequency domain and (complex) Laplace variable domain. In 
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the first case, simple harmonic motion is assumed and an appropriate 
eigenvalue is defined such that a stability boundary results for a realistic 
value of the eigenvalue. The principal advantage of this type of eigen- 
solution is that it can use a wide variety of available unsteady aerodynamic 
theories and/or experimental data which, owing to simplicity, are defined in 
the frequency domain. The principle disadvantage of this type of eigen- 
solution is that it can predict neither actual stability levels nor 
quantitative stabilizing (or destabilizing) trends. A further disadvantage 
of this type of eigensolution is that it provides an aeroelastic description 
of a structural member which cannot be easily integrated with transfer 
function dynamic descriptions of other structural or control elements. 

For the above reasons, the frequency domain eigensolution approach was 
rejected in favor of the more general complex frequency (Laplace variable) 
type. Thus, the basic form of eigensolution problem considered in the G400 
analyses is represented by the following matrix equation: 

w{o}+ [B]{o}+ rc]{o}= {0} 


where the homogeneous solution is given by: 


|q| = [ < i>][exp(\'|')3 


(141) 


(142) 


The mode shapes and complex frequencies are given, respectively, by the 
columns of $ and X_^ (= + ia)^) . The mode shapes are generally complex and, 
like the eigenvalues, occur in complex pairs, when not strictly real. 

Numerical Considerations 


The general solution of Eq. (141) for and X^, as defined in Eq. (142) 
must account for the nonsymmetry of the A, B and C matrices comprising the 
eigenvalue problem. Within the G400 analyses, Eq. (141) is solved using the 
QZ algorithm method described in References 26 and 27 and available as part 
of the EISPACK Matrix Eigensystem Subroutine Package (References 28 and 29) . 
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Principal Assumptions 


The additional eigensolution development undertaken in this study, to 
account for structural and aerodynamic sweep, was accomplsihed within the 
existing mathematical structuring already established for the G400 analyses, 
as reported in Reference 1. Within this context, the following principal 
assumptions were used to guide the additional formulations: 

1. The elastomechanical modeling of the unswept blade, including 

linearization about an initially deflected position, is retained. 

This includes the mode deflection implementation of the AEI nonlinear 

torsion excitation as given by Eq. (30), but only with regard to the 

elastic curvatures due to modal deflection, v M and w "• 

e e 

2* Within the bending equations, the effects of elastomechanical coupling 
due to structural sweep can be formulated using the equivalent prebend 
principle. 

3. Within the torsion equation, the couplings due to structural sweep are 
formulated using the force integration approach of Eq. (40). However, 
consistent with Assumption 2, the loadings required for Eq. (40) are 
evaluated using the approximations for small structural sweep. 

4. The aerodynamic modeling formulated for large aerodynamic sweep in an 
above section is used, as formulated, with appropriate linearization. 

5. The aerodynamic submatrices comprising Eq. (141) can be accurately 
(and more practically) evaluated using direct integrations of the 
various Galerkin weighted partial derivatives (in implicit compact 
form) . This is contrary to the more usual approach of first evaluating 
numerous integration constants and then forming the matrices from 
these constants. 


Prebend Equivalency Principle 

The prebend equivalency principle is a technique whereby existing eigen- 
solution terms formulated for zero elastic axis offset (structural sweeo) can 
be used, with appropriate substitutions, to yield the new structural sweep 
related terms in the bending equations. The basis of the principle derives 
from the definitions for the deflection correction functions, Eqs. (4) and 
(5), discussed in an above section. These deflection correction functions 
are typically formed from integrations of nonlinear combinations of twist 
and bending deflection functions. Hence, the integrands of these integrals 
are of the general form: 


102 


(143a) 


Aw (,) 


where : 



f(x,) = 0'(x,) v e (x,) = (B's + e't) V e (143b) 

The prebend equivalency principle states that the effect of elastic axis 
offset is to augment the integrand as if it were an increment to the elastic 
axis displacement due to bending , but omitting the combination of built-in 
twist with built-in elastic axis offset: 


f + Af = 00 Ve + 0e(ve + Yio ea ) 


(144) 


Perturbations to this equation then become: 


8(f + Af ) = 00 8 vq v + 00 8\% 


+ v e S5 e + yio EA 80e 


(145) 


The interpretations of these equations is twofold* Equation (144) 
indicates that for those (nonlinear) terms, which are proportional to finite 
deflection, analogous terms should arise with elastic axis offset substituted 
for initial bending deflection. Equation (145) indicates that for those 
linear terms involving built-in twist (term (i) ) analogous linear terms should 
arise involving built-in elastic axis offset (term (l\^ ). Note that the 

second and third terms are nonliner and contribute nothing to the equivalency 
principle. Thus, this prinicple must be implemented by examination of 
both the linear and nonlinear terms appearing in the bending equations. 

As per the third principal assumption in this section, the effects in the 
torsion equation are treated in a different manner. 

Nonlinear Source Terms 


In the original eigensolution formulation summarized in Reference 1, 
various nonlinear elastomechanical terms emerged which are of the form: 
qvfc* ( s n and 3w 1 ( S n The quantities q^f and q Wm i , are the initial 

generalized (modal) deflections in edgewise bending and flatwise bending, 
respectively. The general quantities and 6n refer to some appropriate 
integration constant, and to one of the perturbational response variables, 
<5q w ,> ^vfc* (and/or their time derivations), respectively. The prebend 

equivalency principle as outlined above states that built-in structural sweep, 

y 10 , and z^o > should behave as initial modal bending deflections. Thus, 

EA EA 
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in the present context the previously derived elastomechanical terms 

involving q w and q v , can be used in this heuristic manner to formulate 
m 

the first set of terms arising from built-in structural sweep (prebends) . 
Linear Source Terms 


The perturbational results of Eq. (145) pertain to those linear 
elastomechanical terms which are of the form where S£ is some appropriate 

built-in twist (0g) related integration constant and is a bending 
perturbational response (Sq^, 6q w ^, or derivatives). The prebend 
equivalency principle states that the linear terms involving elastic twist 
( 60 e ') perturbation in the presence of built-in elastic axis offset 

( y lO E A^ an< * ^IOea^ should behave the same way as the above S^dg terms. 

The detailed practical implementation of the prebend equivalency principle 
relating to the nonlinear and linear source terms is summarized in Table IV. 
Note that with any entered pair of equivalencies, both the functional and 
its first spanwise derivative are given, were appropriate. 


Perturbational Airload Matrices 

The perturbational airload matrices and submatrices required for the 
Eq. (141) eigensolution format are formulated using principal assumptions 
4 and 5. The actual calculation of these matrices proceeds in various matrix 
multiplication steps representing the steps of a chain-rule differentiation 
procedure. One advantage of this approach is that it maximizes the accuracy 
in the linearization process because it is implicit, compact and, hence, 
reasonably tractable. The major steps in this chain-rule, matrix multipli- 
cation procedure are described in the following subsection. 

Partial Derivative Matrix of Unsteady Airloads 

The appropriate starting point for calculating the perturbational aero- 
dynamic matrices is the collection of formulae defining the total load 
distributions, Eqs. (57) through (61) and (100). This subsection formulates 
the various partial derivative matrices leading to the partial derivative 
matrix of airloads with respect to the intermediary (aerodynamic) perturbation 
vector 6 Z. Specifically, the following partial derivative matrix is sought: 


[DFOZ] = 


_ aq °«, /az i . 


( 146 ) 
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TABLE IV 


SUBSTITUTIONS USED FOR IMPLEMENTING PREBEND EQUIVALENCY PRINCIPLE 


Bending Functionals 

| Structural Sweep Substitution 

Nonlinear 

Source Terms 

Y ^m' y w m ' 

2 7 1 

Z I0 EA i z IO EA 

X V '• V 

^IO EA i ^IOea 

Av ^'l % m '; Av m f q Wm , 

^V EAji A Ve <*' 

Aw ®kj V* AWe kf >/qv k' 

Aw EAj, Aw E if 

K mm ' W Sq Wm + u e Ehk / Qv k ' Sqv k ) 

( Au EA Fm Sq Wm + Au EAE ^Sq Vk ) 

Linear Source Terms 

Au Bm S W Av£,' Bq Wm 

^EAjS^jiAv^'Sqe. 

Aw Bk Sq v k i Aw B k 2)/ SQv k 

Aw eaj s <l 0 j i Aw ea ^ Sq$. 

AV k 8q V|(i Avf'Sq,,, 

AV ea - 8q^. , A/ EA ^Sq^ 

AW m 8q Wmi AW m K|, 8q Wm 

AW EAj 8q 8jl AW E ®8q Sj 
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where p a and Pa zlQ are, respectively, 
edgewise y and flatwise directions: 


the airload distributions in the 




cos0+ p 02a sin© + p 0x8 (^ e cos®+ Z6«sine-AA) 


(147a) 


Pa, = " Po. sin0 + Pa z cos®+ Pa (z' 6e cos0- y' e sin®) (147b) 

Z I0 y 8 “ 


and <5z. are the elements of the intermediary aerodynamic perturbation vector, 

6{Z}. 3 

Pe£turbat_ions[ £f_Basi : c_Al L r loads^ 

The first step in evaluating Eq. (146) is to form the partial 
derivatives of the "8" coordinate system force air loadings. Using Eq. (57) 
through (59) one obtains: 
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Of the three contributions to the perturbational " 8 " coordinate 
system airloads shown in Eq. (148), only the first two perturbational 
vectors have to be further expanded to forms relating to components of 
the 6 Z vector. 

Vector 

The elements of the intermediary perturbation vector fall into two 
main groupings: the first relating to blade kinematics, 6 Z^, and the second 
relating to the Pad£ augmented state variables, 6 Z 2 : 



The twelve (12) elements of Z^ are respectively defined as: 


[z,J = [Up 5 ,UT 5 ,UR 5 >y6 e * z 6e» y6e» z 6e>®Ni^N«y5»Z5» u eJ 


(149b) 


The (maximum of) forty (40) Pade augmented state variables derives from four 
Pade variables defined each for lift and moment (for a subtotal of eight) 
taken each with a spanwise shape distribution given by the first five (5) 
Legendre polynominals . Note that some elements of the Z ]_ vector, as defined 
in Eq. (149) are constant or unused within the scope of the present study 
and therefore have zero perturbational values, effectively. 

Perturbational Section Coefficients 


The next step in evaluating the DFDZ matrix (Eq. (146)) is obtaining the 
partial derivatives of the airfoil section coefficients appearing in Eq. (148) , 
with respect to the element of the intermediary perturbation vector. This 
in turn depends on the type of unsteady airloading selected and the assumptions 
made regarding each: 
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(Quasi-static) : 


c i ) | 

( dcj i /da H 

8c i /8m 

c dp( . < 

1 ^dp/^®N 

ac dp /dM 

Cd ,( ' 0 

acd s /aM 

c m / 1 

( dc m /do N 

dCm / 8M 


fM 

\m n ; 


(150) 


( Pade Approximant) : 

From Eq. (113) the following perturbational form can be written: 

cx 

c d p 

C <*s 
'm 





# 

- 

1 

<p \ 
™N 1 

dc^ / d<t> 

dc z /de N 

0 

8c / 80 N 

dc jL /dl z 

1 

1 e H I 

0 

0 

0 

0 

0 

8 < 

1 mA 

0 

0 

0 

0 

0 


1 

E 

o 

a^m/ ddft 

0 

dCm / d&H 

dc^/d7-2 


M 


(151) 


where 5Z^ again represents the perturbations of the Pade augmented state 
variables, Eqs. (114) through (117). Note that in the case of the quasi- 
static airload option, So^ is equal to the sum of 6<j> N and 60 N - Both 5<p 
and 6Mjyj can be further chain-rule expanded using Eq . (50), (48) and (51): 


S *N = -1 (U T 8U p -Up8U T ) 

Uw 


(152) 


SM n = — (U T 8U t + UpSUp) 


a co 


(153) 


Note also that in the case of quasi-static airloads, the 60^ dependency appears 
in Equation (148) (k ^ 0) , whereas in the case of Pade approximate airloads, 

* (151). 


k = 0 and the 0^ dependency appears instead in Eq 
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Perturbation Section Velocities 


The last major step in evaluating Eq. (146) is to obtain expressions 
for the partial derivatives of the various perturbational section velocities 
appearing in the second contribution to Eq. (148). The appropriate relation- 
ship defining these velocities comes from the above section on aerodynamic 
sweep : 



(154) 


Taking perturbations of this equation leads to the required perturbations 
in U T and U^. The remaining perturbations in and U come from 
differentiations of Eqs. (48) and (41): 

SUn= i [u T SU T + UpSUp] ( 155 ) 

8U = JJ [u t 8U t + U p 8Up+ U r 8Ur] (156) 


Perturbational Equations for Pade Augmented State Variables 


Using Eq. (119) as a basis, perturbational forms 
equations for the (lift) Pade augmented state vectors 

of the differential 
can be written as: 

8x -(^)8x = 8^ 

(157a) 

Sy - '(•£■ fr 2 )Sy = 84> 

(157b) 

8z-(^r Pj)8z = 8l N 

(157c) 

Sw-(-£pJ8w = 80 N 

(157d) 
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These four differential equations, identical in form, for either the lift or 
moment perturbational augmented state variables strictly apply at only one 
spanwise station. For the (NSEG) spanwise stations, there would, for complete 
accuracy, have to be (8 x NSEG) differential equations (4 for lift and 4 for 
moment) to be solved simultaneously with the blade modal equations to complete 
the eigensolution formulation. For this most general case, it appears that a 
large eigenvalue problem would have to be solved. A practical method for 
keeping this problem tractable is discussed in the next subsection. 

Reduction ^f_Ajagmeiitjed_S_ta_te_yarl L ab^le^ ^ej^r^e^-o^fj^F_ree^dom 

The (8 x NSEG) degrees-of-freedom due to distinct augmented state 
variables at each spanwise station can be reduced by expanding the augmented 
state variables in terms of Legendre polynominals defined along span r. For 
example: 


N 

8x(r,<M= I Pj _ , (r) 8 XjOfO 


(158) 


where the Legendre polynominals are defined as: 

P 0 (»’)= I P 3 (r) = 20r3 - 30 r 2 + I2r - I 

FJ(r) = 2r — I (^(r) = 70r 4 -i40r 3 + 90r 2 - 20r + I 

P 2 (r) = 6r 2 — 6r + 1 

Substitution of Eq. (158) into the differential equation for 6x, Eq. (157a) 
yields: 


* N 

Sx. - (2i -l) I <rJj X> 8Xj*(2i-l)ftf,(*) 

(159) 

jC ( ^ (r)P i-i (r)P j-. (r)dr 

(160a) 

8 f jW= r ,p i ( r )8<?(r»dr 

(160b) 
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Recalling the definitions for the inflow angle, <J>: 


<t > = tan 1 (Up / Ur) 


and the normal flow velocity, U : 

N 

2 2 2 

U N = Up + U T 

one can easily shown that Eq. (159) becomes: 




+ — (UjUp 2 - UtUj 2 — 2UpUpUx)8Up 

Uu 


(161) 


Un 


, (Up Up 2 - UpU T 2 + 2UTU p u T )8UTjdr 


Similarly, one can show the following: 


8y - (2i-i)£ 0 jy) £y. = [same RKS os eq. I6l] 

j=i iJ 1 *• J 


(162) 


where : 


C7. 


<y> 


•m 


z ( '» p i 


-i ( r) p J-|(rJdr 


(163) 


Substituting similar Legendre polynominal expansions into the 6z and 6w 
equations in Eqs. (157c) and (157d) , and recalling the modal expansion of 
pitch angle perturbation: 

NTM 

S0 N = Z Tg. Sqn 

j J J 
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one can easily show the following: 


SZ. - (2i-|) 


N 

i 


j=« 


V ™ 82. = 


(2 i -l)IJo' drP i-i lr) %i (r)8 ^i 
j J J 


(164) 




same RHS as eq.164 




(165) 


where: 

ov(z) = J' (-Jftyr) Pj_ ( (r) Pj_,(r)dr (166a) 

r Jo (cl p 2 lr,P i-l (r)P i-' (r,dr <166b) 


The smaller set of first order differential equations, Eq. (161), (162), 

(164) and (165) , for optional choice of N in the Legendre polynominal 
expansion, Eq. (158) can be used to replace the larger the set of differential 
equations defined by Eq. (157). 


Decomposition of Intermediary Perturbation Vector 


All of the above development has been directed to the evaluation of 
Eq. (146) and to the formulation of the Pade augmented state variable 
equations, both of which are defined in terms of the indermediary perturbation 
vector, 6{Z). However, the ultimate eigenproblem formulation must be in 
terms of perturbations of the blade generalized (modal) coordinates (<5q w , 

6Zo. The method for 


6q, 


>V k and 6qg^ ) and the Pade augmented state variables, 
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eliminating the 6{Z} vector and introducing in its place these generalized 
coordinates is an appropriate second major chain rule, matrix multiplication: 


S{ z H t ,]{SQ„} + [t 2 ]{SQ^ ♦ [t 5 ]{So4 


( 167 ) 


where the blade coordinates are given by: 


UqJ * sUw,... q WNFM ... q„... q» NEM ...q 9| ...q eNIM J (168a > 


and where the weighted Pade augmented variables are given by: 


(168b) 


The details of the T^, and T partial derivative matrices are obtained 
using the formulation given in the above section on structural twist and 
sweep, and the above subsection relating to the perturbational forms for the 
Pade augmented stall variables. While these formulations involve straight- 
forward partial differentiation, they are quite tedious and presentation of 
the details herein would contribute little to the clarity or usefulness of 
this report. 

General Application of the Galerkin Method 

The remaining formulations required to define the final matrix operations 
for calculating the required aerodynamic matrices are the approriate Galerkin 
weighted integrations. Define the generalized (modal) excitation vector 
commensurate with the vector of generalized coordinates: 



where, by Eqs. (146) and (167): 
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( 170 ) 


8 tel = 

and where the Galerkin weighting matrix is given by: 



(rWj7 AWj) 

-Aw k 



A Yi 


(rv k -;Av k ) 



0 

0 


A 



(171) 


Finally, the three required perturbational airload matrices can be written 
as : 


[bAI|]=— / 0 '[r B ][DFDz][r,]dr 

(172) 

[cah]= -/ o '[r B ][ DFDZ ][ T2 ] dr 

(173) 

M = - /„' [r B J[DFD z ][T 3 ]d r 

(174) 

For completeness, the following additional matrices 
Pade airloads modeling can be defined: 

appropriate to 

r i r **i 


|AA2lJ= [df/dQ B \ 

(175) 

[ba2i] = [df/dOe] 

(176) 

M*M 

(177) 


where the f and a vectors are ensembles of similar constants defined as by 
Eqs. (160), (163), (164) and (166). 
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Structuring Techniques for Semi-Canonical Form 


The purpose of this subsection is to integrate the definitions for the 
various components submatrices given above into the semi-canonical form 
required by the EISPACK matrix eigensolution algorithms. Specifically, 
the general form of the perturbational equations of motion, Eq. (141), can 
be written in its most detailed form: 


> n 

o 

1 

JSL 

Bll + BAII 1 0 _ 

1 


r 

AI2 0 

W 

BA2I I 



+ 

CII+ CAII 

CAI2 
^ 

8 \ 

kl 


0 

CA22 

1 

l Qp J 



(178) 


The elastomechanic submatrices. All, Bll and Cll are obtained from Reference 1 
and as outlined above regarding the prebend equivalency principle. The 
remaining aerodynamic matrices, BA11, CA11, CA12, AA21, BA21, CA22 are 
given by Eqs. (172) through (177). 


The three semi-canonical forms which are now required are those for 
the vacuum, nonvacuum- quasi-static, and nonvacuum-Pad4 airload optional 
calculations . 


Vacuum Case 


For this case, all aerodynamic matrices are deleted and the resulting 
semi-canonical form becomes: 


“ 





All 1 Bll 


0 1- Cl 1 



— ^ 

X “ 

\ — — 



0 | I 


O 

H 



L 


L ^ 




(179) 


Nonvacuum, Quasi-Static Airloads Case 


For this case, only the principal upper left-hand corner aerodynamic 
submatrices are retained and the semi-canonical form becomes: 


All | Bll + BAII 

o" r i 


X- 


0 I - Cll - CAII 

T| 3 



(180) 
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Nonvacuum, Pade Approximant Airloads Case 


For this most general case, the semi-canonical form becomes: 








All 

0 

| BII+BAII | 0 

r i jo 

x- 

0 | —CM - CAM | -CAI2 

Tj 0 0 

8< 

1 U ® 1 

AA2I 

I BA 21 |T 

■a 


0 | 0 | CA22 

1 



( 181 ) 


Note that the maximum dimension of is 10 and that of Q is 40; 
thus, the total maximum size of the matrix eigenproblem defined by Eq. (181) 
is 60 x 60. 
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SAMPLE CALCULATIONS MADE USING G400PROP 


The SR-2 Model prop-fan blade was chosen to provide a realistic physical 
data base for the purpose of demonstrating the capabilities of the G400PROP 
analysis. This selection was made for two reasons: First, since the SR-2 

blade is unswept, complete beam data exists for this prop-fan configuration 
(in contrast to other, swept prop-fan configurations for which accurate 
equivalent beam data do not yet exist). Second, results from NASA conducted 
stall flutter tests of this blade exist and were readily accessible. 

The appropriate SR-2 blade beam data were assembled and input to the 
G400PROP analysis, first to determine blade natural frequencies (E159 
preprocessor) and then to evaluate coupled mode characteristics and 
aeroelastic stability. The following subsections describe in turn the blade 
configuration and selected operating conditions, the uncoupled modal 
frequency calculations, blade coupled frequencies and mode shapes as 
calculated by the eigensolution (both without and with sweep), and finally 
correlation results comparing experimental stall flutter characteristics with 
the G400PROP predictions. 

Description of Selected Blade Configuration 
and Operating Conditions 

The SR-2 prop-fan propeller model is of solid steel construction, has 
a .6223 meter diameter, and is configured with eight "shovel tipped" blades 
(no sweep). The planform of this model design is shown in Figure 15, and 
a summary of the pertinent geometric and other measured parameters is given 
in Table V, Also included in this table are the various dynamic parameters 
which were either calculated or estimated. The blade uncoupled modal 
frequencies listed were obtained using the E159 preprocessor portion of the 
code and the frequencies are presented in units of both Hz and per rotor 
frequency (P) , based on the given design tip speed. Rough estimates of 
the viscous equivalent structural damping values were estimated on the 
basis of the stall flutter results; these estimates are discussed in 
greater detail in a subsequent subsection. The torsion stress/torque 
coefficient, x /T, at the 19.05 cm spanwise location was calculated using 
blade geometry and appropriate formulas from Reference 30. The nominal 
section properties for the SR-2 model blade are then given as functions of 
nondimens ional spanwise location in Figures 16 through 24. 
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Figure IS. Planform of SR-2 Model Prop-Fan Blade 
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TABLE V 


SR-2 MODEL PROP-FAN PHYSICAL PARAMETERS 


Design Parameters 


Model Values 


Tip Speed, ftR, m/s 

Rotor Speed, ft, rpm 

No. of Blades, b 

Radius, R, m 

Solidity, o 

Blade Root Offset, e 

Preconing, 6 , deg 
B 

Prelead-lag, 6 deg 
B 

Fabrication Material 
Parameters Calculated or Estimated 
LJncou.p_l e d_ Mo_d e _N^t ur a_l_F r_e^u en£i_e s 


277.01 

8500 

8 

.3112 
.565 
. 1567R 
0.1 


4340 stainless 
steel 


1st Flatwise Natural Frequency, u) , Hz 


233.34 (1.577P) 


2nd Flatwise Natural Frequency, u) , Hz 

W 2 

3rd Flatwise Natural Frequency, u) , Hz 

W 3 

1st Edgewise Natural Frequency, 0) , Hz 

V 1 

1st Torsional Natural Frequency, u> , Hz 


541.94 (3.825P) 


1037.74 (7.325P) 


1030.59 (7.274P) 


627.64 (4.430P) 


2nd Torsional Natural Frequency, u> , Hz 


1246.79 (8.801P) 


^t£u£tural_C£i_t ic_a_l Damping Ratios 


Flatwise Modes 
Edgewise Mode 
Torsion Modes 


0.008 

0.008 

0.008 


Torsion S^tr^e^s^P ^t^hin^ Moment 


t /M (@ r «= 19.05 cm) 


(63 . 77/in ) 


120 


BLADE OFFSET 



LO 

CM 

O 


o in o ^ 

CM r- n- o 

• • • . 

O O O o 

SNOiinaimsia H3i3wvwd nouoss 3avig 


o 

o 


o 

d 


121 


NONDIMENSIONAL SPANWISE STATION, 
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NON DIMENSIONAL SPANWISE STATION, F 

Figure 17. Weight and Inertia Distributions for SR-2 Blade 
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Figure 23. Twist and Fictitious Elastic Offset Distributions for SR-2 Blade 
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Figure 24. SR-2 Blade Section Twist Stiffness Distributions 



The seven cases selected to define operating conditions for the SR-2 
prop-fan are described in Table VI. The first two conditions are defined 
by a blade pitch angle and rotor speed representative of unstalled operation, 
both without and with a fictitious (but representative) structural sweep 
distribution. These two cases were selected solely as vehicles for 
demonstrating the capability of the code and are not intended to show 
any quantitative performance characteristics of the SR-2 blade. While 
case 1 could potentially be correlated with test data, if such data were 
available, case 2, with the fictitious structural sweep of 35.3 deg, can 
only have meaning in a theoretical sense. The purposes of these two cases 
are, therefore, to demonstrate program operation, as check cases, and to 
establish qualitative sweep related trends. 

The fictitious sweep selected for case 2 was arbitrarily constructed 
to have a sweep angle of approximately 35.3 deg at and outboard of the 75% 
span location with parabolic variation inboard to the offset location. This 
fictitious sweep (elastic axis offset) distribution is shown in Figure 23. 

The forward flight speed was chosen to maintain an unstalled blade angle-of 
attack of about 4° at the .75R spanwise location. The number of flatwise, 
edgewise, and torsion modes used were chosen to exercise the eigensolution 
to its fullest extent. 

The last five conditions shown are identified for the nominal (nonswept) 
SR-2 blade and are intended for correlation with -the stall flutter test 
results, as discussed in a later subsection. With the seventh correlation 
case, the edgewise mode was deleted due to the expected inactivity of this 
relatively high frequency mode, and only the first torsion mode was retained 
commensurate with the predominantly first torsion mode excitation test results. 

The variations in uncoupled blade modal frequencies with tip speed are 
presented in Figure 25. These frequencies, of course, lack the coupling 
effects of twist, precone, prelead, etc., which the G400PROP analysis provides. 
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SELECTED OPERATING CONDITIONS FOR SR-2 PROP-FAN 
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BLADE NATURAL FREQUENCY, 



ROTATIONAL SPEED, RPM 


Figure 25. Variation of SR-2 Blade Natural Frequencies with Tip Speed, Uncoupled Modes 


132 


PER REV LINES 



Results of Eigensolution 


Both cases 1 and 2 were run using the eigensolution analysis to establish 
satisfactory operation of the code in predicting realistic coupled mode 
vibration characteristics as well as reasonable aeroelastic stability estimates. 
The results of these eigensolutions are presented in complementary forms in 
Table VII and Figures 26a, 26b, and 26c. The table summarizes the eigenvalue 
results for both the unswept and swept configuration. These consist of the 
(vacuum) coupled made frequencies, as well as the (nonvacuum) aeroelastic 
stability damping levels. The three portions of Figure 26 summarize the 
vacuum case eigenvector results for these two configurations for the six 
coupled modes. The results are presented in phase-plane format for the three 
components of tip motion: inplane motion (y^), out of plane motion (z^) and 

pitching (0). For each of these modes the results for the nominal (unswept) 

SR-2 are compared with the corresponding results for the fictitiously swept 
blade. 

Figure 26a shows comparisons of modes 1 and 2 for the swept and unswept 
configurations. In particular, it can be seen from the phase-plane diagrams 
that out-of-plane motion is the dominant factor in modes 1 and 2, with 
significant pitching motion only in mode 2, for both the swept and unswept 
configurations. This is to be expected since in both configurations 
mode 1 represents primarily the first flatwise mode (IF), and mode 2 represents 
the second flatwise mode (2F) with some first torsion mode (IT) coupling. 

The opposite signs of the inplane (Y5) and out-of-plane ( Z 5 ) phase-plane 
vectors indicate that the bending is primarily flatwise. Also, since the 
swept blade second mode pitching vector is relatively in phase with the 
bending, it can be seen that sweep is generating a strong coupling between 
torsion and flatwise bending as is expected. 

Figure 26b shows similar comparisons for modes 3 and 4. The phase- 
plane vectors illustrate a great deal of pitching motion for modes 3 (IT) 
and 4 (3F, IE, 2T) , but a significant amount of inplane and out-of-plane motion 
for both blade configurations only for mode 4. The bending in mode 4 is 
similarly indicated to be primarily flatwise. Most 3 has a coupled frequency 
which closely resembles that for pure first torsion mode (IT); this mode is 
essentially the coupled form of the first torsion mode. Mode 4 has a frequency 
which is very close to that for third flatwise (3F), first edgewise (IE), and 
the second torsion (2T) modes. 

Figure 26c shows comparisons of mode and 6 for the two blade 
configurations. Because the Y5 and Z5 components of mode 5 are now in-phase, 
the bending is indicated to be edgewise. The 5th mode for the swept configu- 
ration unexpectedly exhibits a large amount of coupled pitching motion with 
the second torsion mode (2T) . Mode 6 is seen to represent pure second torsion 
mode (2T), and is thus the coupled form of this mode. 
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eigenvalue: X « a ± iu> (nondimensional with respect to Q) 
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Unswept and Swept Blade Configurations, Modes 1 and 2 
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Figure 26b. Coupled Mode Shape Comparisons, at the Tip Section, for 
Unswept and Swept Blade Configurations, Modes 3 and 4 
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Stall Flutter Correlation Cases 


Figure 27 summarizes the experimental and analytic stall flutter 
results for the statically thrusting SR-2 Prop-Fan model blade. The 
experimental results are presented in the form of isostress curves for the 
1/2 peak-to-peak (1/2PTP) torsion stresses measured at the 19.05 cm spanwise 
location, as shown in Figure 15. These results were obtained from stall 
flutter tests on model Prop-Fan configurations.* The reduction of the test 
data to 1/2PTP torsion stresses was accomplished using the manufacturer's 
quoted gage factors for the torsion strain gages used. 

The analytic G400PROP calculations consist of cases 3 thru 7 given in 
Table VII and are indicated on Figure 27 by the open or closed square 
symbols. These five cases were selected with blade pitch angle-rotor speed 
combinations which would appear to result in both stable and unstable 
calculated responses. Cases 3 and 7 were selected to overlap partially 
with case 5 in either blade pitch angle (Case 3) or rotor speed (Case 7). 
Case 5 was deemed to be a representative condition strongly associated with 
stall flutter, as opposed to the deep stall buffeting which was observed at 
significantly higher pitch angles. Correspondingly, Cases 3 and 7 were 
deemed to be probably stable conditions. Furthermore, in the course of 
the study, additional conditions. Cases 4 and 6, were introduced to enhance 
the definition of the apparent stall flutter boundary. In all the G400PROP 
calculations, the dynamic stall parameters appropriate to the NLR-1 airfoil 
were used. 

In the course of performing the G400PROP calculations, it was noted 
that the degree of stall flutter response obtainable was a strong function 
of the values selected for structural damping. With a sufficiently high 
value of damping (0.02) the blade was stable and a stall flutter condition 
could not be induced. The final value of damping used (0.008) was selected 
on the basis that it produced "stable" limit cycle oscillations which 
neither grew nor diminished once the instability "locked— in.' 

The final G400PROP calculation results for the five conditions are 
shown in Figure 27 either with an open or closed symbol denoting stability 
or instability, respectively. Where appropriate, the calculated 1/2PTP 
torsion limit cycle torsion stresses are shown parenthetically. Generally, 
the detailed results of this figure show the correlation between experiment 


* These tests were performed under contract NA53— 22755 and are summarized in 
UTRC Report R81-335414, "Static Stall Flutter Tests of HSD Prop-Fan Models." 
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Figure 27. Stall Flutter Correlation Results for the SR-2 Blade and G400PROP Calculations 



and analysis to be quite good in that: 

(1) Stall flutter was predicted for conditions close to those where 
it was measured, 

(2) Stable conditions were predicted where the propeller was found 
to be stable, and 

(3) The maximum 1/2PTP stress levels correlate well with those 
measured (despite the disparity in methods used for obtaining 
the stresses, experiment vs. analysis). 

It would appear that the principal disagreement between the experiment and 
theory is the pitch angle at which the apparent stall flutter boundary 
occurs; a disparity of about 6 degrees is indicated. Possible reasons for 
this disparity are as follows: 

(1) The inadequacy of using uniform inflow for the statically 
thrusting condition. Use of the PANPER code to obtain a more 
realistic flow field for the static thrust condition is 
practical and could have been used, but was outside the scope of 
the present study. 

(2) The uncertainty in the static stall characteristics of the NACA 
16-series airfoil data used. 

(3) The unknown impact of cascading effects on the airfoil stall 
characteristics. 
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GENERAL DESCRIPTION OF THE G400PROP PROGRAM STRUCTURE 


The G400PROP Rotor Aeroelastic Analysis is an extended, special purpose 
version of the parent G400 analysis described in Reference 1. For applica- 
tion to the special requirements of propellers, the G400 analysis was ex- 
tended to include an internal preprocessor for uncoupled beam normal modes, 
the elastomechanics accruing from structural sweep, and an upgraded, more 
rigorous aerodynamic modeling. This aerodynamic modeling includes aero- 
dynamic sweep, stalled and unstalled unsteady effects, and three-dimen- 
sional variable inflow. This version now includes an expanded eigensolution 
and direct coupling with the PANPER variable inflow analysis. The program 
was coded in FORTRAN IV and developed on the UNIVAC 1110/81A computer system. 


Program Structure 

The G400PROP program is structured in a generally conventional hierarch- 
ical fashion with an initial routine (MAIN) , appropriate modular elements 
and a collection of utility subroutines and/or functions. Figures 28 through 
31 present flowcharts showing the four sequential portions of the G400PROP 
program structure; each of these figures highlights one of the four major 
portions of the code. Note that the (A) end terminal depicted in Fig. 38 
corresponds to the beginning terminal of Fig. 39, and so on with terminals 
( 5 ) and (c) with Figs. 29 through 31. In these figures, where appropriate, 
the names of the subroutines performing the functionally labeled operations, 
are given parenthetically above their respective functional symbol. 

The hierarchical subroutine calling structure is presented in tabular 
form in Table VIII; this functional structuring also serves to list the 
elements (subroutine, functions and/or data), contained in the G400PR0P 
code. Note that the hierarchical structure is alternately indicated using 
indentations with symbols denoting level and/or with parentheses. Note 
also that this table lists the routines in alphabetical order within any 
given hierarchy and not with respect to calling order, size or importance. 
Following this table is Table IX which presents a common block/subroutine 
cross reference list. Finally, after these two tables the last subsection 
presents a brief description of each of the G400PR0P elements. 
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Figure 28. Flowchart for G400PROP — Input and El 59 Preprocessor Portions 
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Figure 29. Flowchart for G400PROP-Part I, Eigensolution Portions 
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Figure 30. Flowchart for G400PROP-Part II, Time-History Solution Portions 
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Figure 31. Flowchart for G400PROP Part III, Transient Spectral Stability Analysis and 
Program Termination Portions 
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TABLE VIII 


HIERARCHICAL LISTING OF G400PROP ELEMENTS 


• MAIN 


■ E159 


ABMFIT 
▲ ELAST 
A INVERT 
ANFMS 


A BMEVAL 
A BMFIT 
4 EIGENE 


▼ INVERT 

▼ SIMUL 


A PCHDAT 
A PCHM0D 
Atmss 

■ INITT 

■ L0ADER 

■ M0DEIN 

■ M0TION 

A ENDC0N 


4 PNPRSU 


▼ BLIN5 
▼DEFLEX 

▼RSPNSS (QQPSET) 

AQSTHRM (HARM) 

Aresetq (SETVAL) 

4 SETVAL 


146 


TABLE VIII (Coat'd) 


A INISH 

4admc0f 

^DPCHEK 

Amajitr 

4aerprf 

^bladel 

▼rspnss (QQPSET) 
▼ setup 
▼spans 


^ALFD0T 

^ALWC0M 

^BLIN5 

: DEFLEX 
HYSDMP 
+ PHID0T 

♦ REVPAD (PADC0F) 

^SHLDM 

♦UNSTCF 


| C0EFF3 
I SYNTH3 

▼ BLIN5 

▼ GETCDS (INTERP) 

▼ GETCLS (INTERP) 
V GETCMS (INTERP) 
f SHLDM 


4CR0UTS 
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TABLE VIII (Cont'd) 


APANM0M (MOMNTM) 

Atable 

^TMARCH 

▼ ADAMS 

▼ EXTRAP 

▼ QPPCAL 

▼ QPPTST 

▼ STRSSS 


- A UNDATA 

■niam 

AAER0 

A DEFLEX 
f LNAER0 


▼ BLIN5 

▼ SHLDM 


MATRXT 

RSPNSS (QQPSET) 


A AIRPAD 


4 DEFLEX 
^FLGDR 

Alnaero 


▼BLIN5 

▼SHLDM 


MATRXT 

PADC0F 
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TABLE VIII (Cont'd) 


^RSPNESS (QQPSET) 

▲ DFYZEA 
▲ EIGSLN 


FILL2 

S0LVE 


▼ EIGRMC 


^RGG 


Iqzhes 

Iqzit 

|QZVAL 

Iqzvec 


Tphase 

▼VECTRS 


^DEFLEX 

^GJR 


▲FILL 

^RSPNSS (QQPSET) 

▲INTC0F 

▲SPNWIZ 

▲SPRINT 

▲TC0UPL 

▲TW0F 


■ PRNT 
■REDATA 

■resetg 

■RESETQ 

▲ SETVAL 
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TABLE VIII (Cont'd) 


■ VARINF 
■WAKUCZ 

Adiscrt 

Afftgen 

Alfct (nprm) 

AREVERS (SBSCRP) 

Asctab 

$SUBS (SBSCRP) 

Alinfit 

Amaximz (IDSCRT) 

A M0DULS 
A SEARCH 
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Element Common Block Names 
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Part 2. Element Name Per Common Block 
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PEESP; C0EFF3 SYNTH 3 

PL0VAR: WAKUCZ 

PIAJS: C0EFF3 SYNTH3 

PREDM: C0EFF3 UNDATA 
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Brief Description of Program Elements 


ADAMS 


This function implements an algorithm for time integration of any of 
the time dependent variables. 

ADMCOF 

This subroutine evaluates the coefficients required for the time 
integration algorithm. 

AERO 


This subroutine generates the aerodynamic damping and stiffness ma- 
trices for the eigensolution assuming a quasi-static formulation, blade mo- 
tion and aerodynamic sweep effects. 

AERPRF 

This subroutine completes the calculations for and outputs the summary 
of aerodynamic performance quantities. 

AIRPAD 

This subroutine generates the aerodynamic inertia, damping and stiff- 
ness submatrices for the 'eigensolution assuming a Pade formulation, blade 
motion and aerodynamic sweep effects. 

ALFDOT 

This subroutine calculates the aerodynamic A parameter using backward 
differencing on the inflow angle and direct knowledge of the time derivative 
of pitch. 

ALWCOM 

This subroutine calculates the unsteady decay parameter, a^, required 
for the unsteady stalled and generalized Wagner function airloads method- 
ologies . 

BLADEL 

This subroutine provides the computational loop structuring over number 
of blades in forming the blade response equations. This loop for G400PROP 
is presently degenerate and the upper limit on this loop is 1. 
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BLIN5 


This subroutine does a tri-variant table look-up of the airfoil section 
coefficients. The three variables are angle-of-attack, Mach number and span- 
wise station. 

BMEVAL 

This subroutine evaluates the bending mode shape and its derivatives 
at spanwise locations other than where they are calculated in the E159 
eigensolution. 

BMFIT 


This subroutine performs a functional (polynomial) fit of the bending 
mode shape for use in subroutine BMEVAL for evaluating bending mode shapes 
at nonstandard spanwise locations. 

C0EFF3 

This subroutine calculates various coefficients needed for the UTRC 
stalled unsteady airloads methodology. 

GROUTS 

This subroutine is a compact simultaneous equations solver used for 
nonteetered rotor configurations. It uses the Grout Reduction method des- 
cribed in Reference 22. 

DEFLEX 

This subroutine evaluates the spanwise deflections, slopes, velocities, 
etc. from the modal responses, forms the sweep transformations and, for the 
eigensolutions , forms various deflection partial derivaties. DEFLEX typi- 
cally operates within the spanwise loop of the calling subroutine. 

DFYZEA 

This subroutine forms the spanwise derivative of the structural sweep 
angles form the input sweep changes per segment length, or from numerical 
differentiation. 
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DISCRT 


The purpose of this subroutine is to compute the magnitude of the 
Fourier coefficient, CMAG, of a set of time-history data, as part of the 
transient spectral stability analysis (TSSA) . 

DPCHEK 

This subroutine tests the input integration steps size for acceptable 
accuracy, and automatically decreases it if the value is too large. 

EIGENE 

This subroutine performs the eigensolution of the bending portion of 
the E159 preprocessor for uncoupled blade frequencies and mode shapes. It 
uses the method of determinant iteration. 

EIGRMC 

The purpose of this subroutine is to coordinate the running of the 
standard nonsymmetric matrix eigensolution subroutine, RGG, to the extent 
of organizing the root pairs which are caluclated for the main (G400) 
eigensolution. 

EIGSLN 

This subroutine coordinates the setup of the main (G400) eigensolution. 
Specifically, it controls the calculation of the matrices comprising the 
eigensolution, and the extraction of the eigenvalues. 

ELAST 

The purpose of this subroutine is to calculate the elastic coefficients 
for flatwise and edgewise bending for the E159 eigensolution. 

ENDCON 

This subroutine serves three main functions associated with the com- 
pletion of the Part II time-history solution: (1) it completes the calcula- 

tions for median and 1/2 peak-to-peak stresses, (2) it controls the harmonic 
analyses of responses hub loads and stresses, and (3) it controls the saving 
of end conditions and other data for use in either the PANPER code or sub- 
sequent G400 runs. 
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EXTRAP 


This subroutine effects a "static" solution on any degree-of-f reedom 
whose natural frequency is sufficiently high to approximate the response 
neglecting the twice time differentiated term in that degree-of-f reedom* s 
governing equation. 

E159 


This subroutine controls the preprocessor calculations of the un- 
coupled modal frequencies and mode shapes for flatwise and edgewise beam 
bending, and for torsion responses. 

FFTGEN 

This subroutine is a standard Fast Fourier Transform calculator, and 
is used in the transient spectral stability analysis (TSSA) . 

FILL 


The purpose of this subroutine is to fill in the nonaerodynamic por- 
tions of the mass, damping and stiffness matrices used in the main (G400) 
eigensolution. 

FILL2 


This subroutine arranges the various matrices and submatrices set up 
by subroutines FILL and AERO or AIRPAD into the semi-canonical form required 
for eigenvalue extraction by the standard eigensolution solver, RGG. 

FLGDR 


This function evaluates any of the first five Legendre polynomials. 
GETCDS 

This subroutine provides internally calculated static aerodynamic drag 
data in place of user provided static airfoil data for usage in the un- 
steady stalled airloads calculation. 

GETCLS 

This subroutine provides internally calculated static aerodynamic lift 
data for usage in the unsteady stalled airloads calculation. 


160 



GETCMS 


This subroutine provides internally calculated static aerodynamic 
moment data for usage in the unsteady stalled airloads calculation. 

GJR 


This utility subroutine optionally obtains simultaneous equation solu- 
tions and/or matrix inversions using the Gauss- Jordan Reduction method. 

HARM 


This utility subroutine performs a Fourier (harmonic) Analysis of any 
time-history string of data. This harmonic analysis uses a negative cosine 
and sine definition for the harmonic components. 

HYSDMP 

This subroutine calculates the increment to blade edgewise bending 
moment to account for hysteretic structural damping. This formation of 
structural damping is dependent on edgewise deflection and the signs of 
rate and acceleration, but not their magnitudes. 

INISH 


This subroutine initializes arrays and logic variables, and nondimen- 
sionalizes parameters, as required for the time-history solution. 

INITT 


This subroutine zeros out the various arrays in common "BLOCKS". 

INTCOF 

This subroutine forms the various elastomechanic integration constants 
used in the eigensolution and, to a limited extent, in the time-history 
solution. 

INTERP 

This subroutine is a general purpose linear interpolation calculator. 
INVERT 

This subroutine is a general purpose matrix inversion, determinant cal- 
culator used by the E159 eigensolution. 
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LFCT 


This subroutine finds the prime decomposition of any integer for use 
with the Fast Fourier Transform subroutine, FFTGEN. 

LINFIT 

This subroutine performs a least-square fit to results from the tran- 
sient spectral stability analysis routine, WAKUCZ. 

LNAERO 

This subroutine calculates partial derivatives of airfoil section 
coefficients with respect to angle-of-attack and Mach number. 

LOADER 

The purpose of this subroutine is the loading of the generic loader 
portion of the input data. 

MAIN 


This is the main program element and directs all major portions of the 
solution flow. 

MAJITR 

This routine contains the time marching computation loop used in the 
time-history solution; it calls the specialized calculation subroutines 
needed in this loop. 

MATRXT 

The purpose of this subroutine is the calculation of various partial 
derivative matrices used by subroutines AERO and AIRPAD to calculate the 
aerodynamic portions of the main (G400) eigensolution. 

MAXIMZ 

The purpose of this subroutine is to maximize the magnitude of the 
Fourier coefficient as a function of frequency in the vicinity of an 
identified high response frequency, as part of the transient spectral 
stability analysis. 
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MODEIN 


This utility subroutine inputs the blade bending and torsion mode shapes 
and their derivatives with respect to span. 

MODULS 

The purpose of this utility subroutine is to form the modulus of a 
vector of Fourier transforms* 

MOMNTM 

This function evaluates an empirical function joining the two branches 
of the actuator disk momentum equation across the vortex ring state based 
on a function given in Reference 31. 

MOTION 

This subroutine controls all the major elements of the time-history 
solution. 

NFMS 

This subroutine controls .the calculation of the uncoupled beam bending 
vibration modes within the E159 eigensolution preprocessor. 

NIAM 

This subroutine performs the following functions: 

1. Performs some of the detailed initializations and/or nondimen- 
sionalizations of logic and system parameters. 

2. Calculates the deflection correction function arrays which accrue 
from structural twist and sweep. 

3. Controls the calculation of the vacuum and nonvacuum main eigen- 
solutions. 

NPRM 


This utility subroutine finds the next prime number given the vector of 
previous primes. It is intended for use with subroutine LFACT. 
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PADCOF 


i 

This subroutine evaluates the Pade coefficients using linear 
interpolation, as appropriate. It incorporates and optionally selects, 
the data from either the Jordan data, Davis and Malcolm data, or generalized 
Wagner function sources. 

PANMOM 

This subroutine forms the derivative of momentum induced velocity for 
use in the time-history solution. 

PCHDAT 

This utility subroutine output punches spanwise array elastomechanical 
data from the E159 portion of the program for subsequent optional explicit 
input to the G400 proper part of the program. 

PCHMOD 

This utility subroutine output punches spanwise mode shape data from 
the E159 portion of the program for subsequent optional explicit input to 
the G400 proper portion of the program. 

PHASE 


This subroutine evaluates the phasing matrices used by the eigensolu- 
tion for stability analysis. 

PHIDOT 

This subroutine calculates the time derivative of inflow angles for use 
in subroutine REVPAD. 

PNPRSU 

This subroutine forms the geometric and aerodynamic parameters 
required by the Propeller/Nacelle Variable Inflow Analysis (PANPER) . 

PRNT 


This subroutine provides an echo print output of the Part II input data 
which pertains to the Inertia, Elastic, Geometric and other Operational Data. 
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QPPCAL 


This subroutine calculates the array of response accelerations for output 
print. 

QPPTST 

This subroutine tests the system degrees-of-freedotn for numerical insta- 
bilities. The criteria used to identify such an instability are the occur- 
rence of three sign changes of increasing amplitude in three time steps. 

QQPSET 

This subroutine sets the blade modal responses and modal rates for a 
given blade from the vectors of system response amplitudes and rates, 
respectively. 

QSTHRM 

This subroutine performs harmonic analyses (using subroutine HARM) of 
the blade modal responses, hub shears and moments, and blade stresses, after 
the responses have converged to periodicity. 

QZHES 


This subroutine is the first step of the QZ algorithm for solving 
generalized matrix eigenvalue problems, as required by subroutine RGG. 

QZIT 


This subroutine is the second step of the QZ algorithm for solving 
generalized matrix eigenvalue problems, as required by subroutine RGG. 

QZVAL 


This subroutine is the third step of the QZ algorithm for solving 
generalized matrix eigenvalue problems as required by subroutine RGG. 

QZVEC 


This subroutine is the fourth step of the QZ algorithm for solving 
generalized matrix eigenvalue problems as required by subroutine RGG. 

REDATA 

The prupose of this subroutine is to read the static airfoil data. 
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RESETG 


This subroutine resets torsion mode shapes which might be displaced by 
the formation of the pseudo-torsion mode. 

RESETQ 

This subroutine places selected terminal conditions into an array and 
writes them to a data file for use as initial conditions for subsequent runs. 

REVPAD 

i 

This subroutine calculates the Pade augmented state variables used in 
the time-history solution of the differential equations for unstalled un- 
steady airloads. 

RGG 


This utility subroutine calculates the eigenvalues of the generalized 
nonsymmetrical matrix eigenproblem. 

RSPNSS 

This subroutine performs the following time-dependent calculations: 

1. Forms the blade azimuth angle and various harmonics. 

2. Sets the impressed control angle and its time derivatives. 

3. Sets the modal response variables from various optional sources. 
SBSCRP 

This subroutine finds the mixed radix representation of an integer 
for use in the Fast Fourier Transform. 

SCTAB 


This utility subroutine exponentiates an angle multiplied by the 
imaginary vector, i. 

SPNWIZ 

This function performs a numerical integration between blade section 
centers of a specific integrand type as required for forming the deflec- 
tion corrections functions due to structural twist. 
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SPRINT 


This subroutine outputs (as optionally requested) the spanwise 
integration coefficients. Although most of these coefficients are required 
only for the eigensolution, some are used in the time-history solution. 

STRSSS 

This subroutine calculates the spanwise stresses and integrated hub 
loads optionally using the force integration or mode deflection methods. 

SUBS 


This subroutine computes appropriate subscripts and exponents for the 
Fast Fourier Transform. 

SYNTH3 

This subroutine is a component of the group of elements comprising the 
unsteady stalled airloads modeling wherein the final calculations of un- 
steady lift, drag and moment are made. 

TABLE 


This utility subroutine performs a table look-up and first derivative 
calculation for use in defining the instantaneous control angle. 

TCOUPL 

This subroutine calculates the coupled torsion modes arising from 
optional use of the rigid body degree-of-freedom with the E159 calculated 
elastic (normal) torsion modes. 

TMARCH 

This subroutine controls the solution flow for obtaining the time- 
history solutions. It furthermore tests for numerical instabilities and 
convergence to periodicity. 

TMSS 


This subroutine calculates the uncoupled torsion mode shapes and 
natural frequencies within the E159 eigensolution. 



TWOF 


This function performs a least-square curve fit calculation on blade 
twist to facilitate subroutine NIAM in the numerical differentiation of 
blade twist. 

UNDATA 

This element contains the empirical coefficients required by subroutine 
C0EFF3 for using the unsteady stalled airloads calculation. 

UNSTCF 

This subroutine controls the implementation of the unsteady stalled 
airloads calculation. 

VARINF 

This subroutine serves as the input/output interface with the PANPER 
propeller/nacelle variable inflow program. 

VECTRS 

For any eigenvalue calculated by the main (G400) eigensolution, this 
subroutine calculates the coupled generalized mode shape and the spanwise 
components of the physical mode shape. 

WAKUCZ 

This subroutine performs the transient spectral stability analysis 
(TSSA) for extracting such stability indicators as characteristic exponent 
and time to half amplitude from the time-history solutions (see Reference 
32) . 
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PROGRAM INPUT DESCRIPTION 


The required input to the program consists of the following major 
data blocks in order of loading: 

I. Airfoil Data 

II. Inertia, Elastic, Geometric and other Operational Data 

III. Blade Mode Shape Data 

IV. Components of Variable Inflow 

Details for preparing the data for each of these blocks are given in 
the subsections which follow. An additional subsection provides informa- 
tion relating to efficient management of the data and to the input and out- 
put of restart data using files. 


I. Airfoil Data 

This data block generally consists of one or more sets of tables of 
two-dimensional lift, drag and pitching moment coefficients versus angle- 
of-attack for various Mach numbers for up to five (5) arbitrary spanwise 
locations. Additionally, if unsteady aerodynamics are used, the static 
stall angles and linear coefficient slopes for both lift and pitching 
moment are included in this table. 

Input Format for First Card(s) 

While actual set-up of this data block follows a basic format 
(described below), specific variations are required on the first card(s) 
of this block depending on optional usage. These variations denote whether 
multiple data sets are to be input for respective spanwise locations, a 
single set is to be input for use on all spanwise locations, or an analytic 
representation of the NACA 0012 airfoil is to be used for all spanwise loca- 
tions. Each of these optional usages is described below. 

Mul 1 tJjpj 1 e_S2 > anwl : se_ Air f oi 1 s_ 

For the case of distinct airfoil characteristics being defined at up to 
five (5) multiple spanwise locations, the first card image format is as 
follows: 


card #1A: 


/ 


NA NRCL NRCD NRCM TITLE (optional) 


(412, A72) 



where NZ would be ignored so long as the sum of the absolute values of the 
quantities NRCL, NRCD and NRCM is 4 or greater. The quantities NRCL, NRCD 
and NRCM are, respectively, the number of radial stations for which multi- 
ple C£ , c^ and c m c /4 airfoil data are to be input, as appropriate, each 
with a minimum absolute value of 1 and a maximum absolute value of 5. 
Normally, NRCL, NRCD and NRCM are input as positive integer numbers. The 
program also provides for the optional input of these data as negative 
values, in which case the extensive output point of these data, as part 
of the normal case printout, is suppressed. Note that at least one of 
these three inputs must have an absolute value of 2 or greater. For multi- 
ple spanwise section properties, an additional card, following the one 
described above, is then required, which begins the input of the c^ air- 
foil data: 


card //IB: / NZ(1) RADCL(l) TITLE(optional) (12, F8.0, A70) 


where NZ(1) is the number of Mach number for which groups of c^ data are 
to be read in for the first radial station; RADCL(l) is the nondimens ional 
radial station at which the airfoil data is defined. 

Single Airfoil Description 

For the case of a single airfoil to be used for all spanwise locations, 
a single first card image is input. This card is similar to the card #1A 
described above, except that the quantities NRCL, NRCD and NRCM are input 
as zero (or 


card #1 


the columns are left blank) ; 


NZ (1) 000 


TITLE (optional) 


(412, A72) 


In this case the card is interpreted as the first card of the c^ 
data with the RADCL(l) information omitted (see description above for 
card //IB) . 

Analyst i 1 c_Ad 1 rjEoj L l_Des_crljp_ti£n__ 

For those optional cases wherein the analytic NACA 0012 airfoil option 
is evoked (see "S" array location 63 discussed in the following subsection) 
the first card image must be a single card with blank or zeroed columns 1 
through 8. For this option, the remainder of the airfoil data block of 
data is omitted. 
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Input Format for Subsequent Cards 

For those cases wherein tabulated airfoil data are to be input, the 
subsequent inputs continue the airfoil c^ data initiated with the #1A and/ 
or #1B cards. Thus, the card image set-up defined below is then input 
NRCL times (or only once, if NRCL = 0): 

card //2: 


cards #2+: 


ALSTAL DCLDAO (F9.0, 9F7.0) 

where: J is the number of data entries to be input into each Mach number 
group. N is the number of angle-of-attack -c^ (abscissae - ordinate) pairs 
to be input. Normally, a maximum of thirty-four (34) a-c^ pairs may be 
input; thirty-three (33) pairs are input if the unsteady option is chosen. 

M is the Mach number appropriate to the data group. A(i) are the N angle- 
of-attack abscissae in degrees and CL(i) are the N lift coefficient 
ordinates. ALSTAL and DCLDAO are, respectively, the static stall angle, in 
degrees, and the lift curve slope at zero angle-of-attack, in per degree; 
these items are needed only if the unsteady airloads option, (A)64, is 
invoked with a value of 2. Note that J is fixed-point formatted, but N 
and M are floating point formatted. 

Cards 2 and 2+ are repeated for each successively higher Mach number. 

A maximum of 12 Mach numbers is allowed and the lowest and highest Mach num- 
bers need not define the total working range as the search technique uses 
the boundary data for Mach numbers beyond the range input. Thus, repeated 
data for zero and supersonic Mach numbers are not needed. The lowest Mach 
number input must contain an angle-of-attack range from -180° to 180° or 
from 0° to 180° depending on whether or not unsymmetric airfoil data is 
being input; all higher Mach number data need extend only from -30° to 30° 
or from 0° to 30° in a similar manner. 

The general format described above is repeated for the c^ and c m 
subblocks in that order but with either card image # IB or #1A used, dlpinding 


. .. A(N) CL(N) 
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on whether multiple airfoil section data are or are not being input (and 
used). The static stall angles and aerodynamic coefficient curve slopes at 
zero angle-of-attack are deleted for the c^ subblock. Lastly, note that the 
total storage allocated for the combined airfoil data is 5000 (which is 
somewhat less than 2500 abscissa - ordinate point pairs). The maximum 
storage should be adequate for most applications and can be "budgeted" 
among the various spanwise sections, as appropriate. 


II. Inertia, Elastic, Geometric and Other Operational Data 
Generic Load Blocks 


This data block includes those items used to define the more detailed 
dynamic features and/or those which are most likely to vary from case to 
case. The card image format for these data is as follows: 

^ ZZ NN L DATA(L) DATA(L+1) ... DATA(L+4) (A1 , II, 14, 5F12.0) 


where: ZZ is an alphameric code item having the value "E", "R", "G”, "D", 

"A", "V", "S", or This alphameric input code determines into which 

generic loader block the data in the card image will be input. These 
generic loader blocks store data pertaining respectively to: the ^159 

preprocessor for calculating blade uncoupled frequencies and modes (see 
accompanying descriptive material). Radial distributions, Geometric charac- 
teristics, State Variables (functions of time), and Solution control. NN 
is the number of data items to be input on a given card; NN must not exceed 
5. L is the location or identifying number within the indicated generic 
loader block of the first data item on the card columns 3-6 right adjusted. 
DATA(L+1) represents the various data items on the card image, columns 7-18, 
19-30, 31-42, 43-54, and 55-66, in the floating point format. The loca- 
tions or identifying numbers for the various data and control items within 
each generic loader block are listed below along with definitions and other 
pertinent comments; note that some data locations are intentionally left 
blank. These intentionally blank data locations represent inputs defined 
for other more general versions of the G400 family of aeroelastic analyses 
but are not appropriate to this specialized propeller version of G400. By 
not reusing these data locations, the generic loader management permits 
input compatibility among the various versions of G400. 

E159 Preprocessor for Uncoupled Modes 

Data input to (E) loader block controls the optional internal calcula- 
tion of the required uncoupled blade modal data. This calculation thus 
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serves as a preprocessor for the main body of the G400PROP analysis. Use 
of this preprocessor is optional to the extent that, if available, these 
required modal data and related spanwise distributions can be input explic- 
itly, without recourse to the E159 calculation as is described below. Use 
of the E159 processor within any given G400PROP run, however, enables the 
calculation and direct transfer of these data for subsequent usage within 
the main G400 calculations and would need not be explicitly input. 

Furthermore, in the input descriptions to follow, there are apparent 
redundancies with regard to items in the (E) block and similar ones in the 
(R) , (G) and (D) blocks. The redundancies and possible implied ambiguities 
where they occur are resolved within the program by a deference to the 
appropriate quantities which are output from the E159 preprocessor. Thus, 
irregardless of what is initially loaded in the (R) , (G) and (D) generic 
loader data blocks, those quantities output from the E159 preprocessor (and 
arising from the input (E) loader data) are the ones which are used and sub- 
sequently ouptut in the echo listing of the input loader data. Note that 
in the descriptions to follow, those quantities which are redundantly defined 
in the (R) , (G) and (D) loader blocks are dagger tagged (t) . In the event 
that use of the E159 preprocessor is erronously requested together with 
explicit input of the modal shape data (see Section III below) the program 
will detect the ambiguity, print a warning and stop. 
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E159 Preprocessor Data 


(E) Location 

Item 

Description 

i 

R 

Rotor radius, ft 

2 

e 

Offset distance of start of elastically 
deformable portion of rotor blade, 
nondimens ional with respect to R 

3 

NSEG 

Number of segments into which the 
blade is partitioned 

4-23 

A r 

Blade breakup distribution of segment 
lengths, in order from root to tip, 
maximum of 20, in. 

24 

NFMC 

Number of flatwise bending modes to 
be calculated 

25 

NEMC 

Number of edgewise bending modes to 
be calculated 

26 

NTMC 

Number of torsion modes to be calculated 

27 

(Control) 

Make nonzero (1.) to output punch loader 
data and uncoupled mode shape data for 
direct input in subsequent runs 

28-47 

(I/c) F 

Section moduli for flatwise bending 
for sections defined in locations 
(E) 4-23, root to tip, in . 

48-67 

(I/ C ) E 

Section moduli for edgewise bending 
for sections defined in locations (E) 
4-23, root to tip, in . 

69 

n r 

Rotor tip speed, ft/sec. 

70 

(Control) 

Make nonzero (l.)to output intermediary 
E159 calculations for debug purposes. 


i 
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E159 Preprocessor Data - (Cont'd) 


(E) Location 


Item 


Description 


71 


100 


101-150 


252 


Au> 


FE 


N 

w 


(w\ As) 


N 


I 


F 


Frequency scan interval for flatwise 
and edgewise bending frequency 
iteration; default value is .5, 
rad/ sec . 

Number of elements in weight 
distribution, table (equal to twice 
the number of constant value spanwise 
increments used). 

Weight distribution table, taken as 

sequential ordinate-abscissa pairs, 

root to tip. Generally, this and 

all other distribution tables should 

include an extended final (tip) 

distribution spanwise increment so 

that the sum of the spanwise increments 

and the offset exceeds the radius by a 

sufficient amount to preclude diminution 

due to numerical round-off. U (w 1 ) =* 

lbs/ in; U ( A s ) = in. * 

k 

Number of elements in flatwise area 
moment of inertia Table. 


253-452 


453 


(Ip, As) 


N 


I 


E 


Flatwise area moment of inertia table 
(sec above remarks for w 1 , As table); 
u ( I ) “ in\ U Us ) = i n , 

F k R 

Number of elements in edgewise area 
moment of inertia table. 


454-653 


( V is) 


655-674 


(x. ) 
l NS 


Edgewise area moment of intetia table 
(see above remarks for w 1 , As table); 
u (I ) = in\ U ( A s ) ■ in. 

E k k 

Stations other than those defined by 
the input breakup ((E)4-23) for 
purposes of evaluating the mode shapes. 


175 



E159 Preprocessor Data - (Cont'd) 


(E) Location 
676-695 

696-715 


716-735 

736 

737-886 

887 

888-1037 

1038 


Item 

m. 

i 



K 


E. 


l 


N 


GJ 


(GJ, A s) 



(k , £ s) 


0 


Description 

Discrete incremental masses added to 
centers of blade segments (defined 
by inputs ((E)4-23), Ib-sec^/ft. 

Explicit flatwise hinge spring rates 
at inboard ends of selected blade 
segments to replace values obtained 
from input flatwise bending stiffness 
considerations, lb-ft/rad. (Note: (E) 
696 and (E)716 should both be 
sufficiently large values to 
approximate infinite stiffness 
retention of a cantilever beam root 
boundary conditions . ) 

Explicit edgewise hinge spring rates 
at inboard ends of selected blade 
segments, to replace values obtained 
from input edgewise bending stiffness 
considerations, lb ft/rad. (See note 
above) 

Number of elements in the torsion 
St. Venant torsion stiffness table. 

Torsion stiffness table (see above 
remarks for w* , As table); U (GJ ) = 
lb-in^; U(As ) = in. 

Number of elements in spar area radius 
of inertia table. 

Spar area radius of gyration table 
(See remarks above for w 1 , As table); 

U (k ) = in. , U( A s ) = in. 

\ k 

Number of elements in torsional inertia 
distribution table. 
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E159 Preprocessor Data - (Cont'd) 


(E) Location 
1039-1188 

1189 

1190 


1191 

1192 

1193 

1194 


Item 

( V 4s) 


(Control) 



( Control) 


A w 

T 


E 

BSCALE 


Description 

Torsional inertia distribution table 

(See remarks above for w*, As table); 

U (I ) = lb-sec ). U ( As ) = in. 

0 k 

Default value implies cantilever 
attachment at the root. Make (-1.) 
to activate root torsion spring 
restraint (input item (E) 1190). 

Torsion root spring rate, lb-in/rad. 
(Note that this input and (D) 34 define 
a redundant capability. This input 
item includes the effect of root 
flexibility in the torsion mode 
directly, whereas (D) 34 is used 
to approximate this coupling 
"after the fact". Consequently, 
input items (E) 1190 and d(34) 
should not both be nonzero. 

Innermost segment which is elastically 
active in torsion; default is 1. 

Frequency scan interval for torsion 
frequency iteration; default value 
in 3 . , rad/ sec . 

Modulus of elasticity; default value 
is 1 x 10 7 , lb/in 2 . 

Factor used to scale the stiffness 
matrix to avoid erroneous zero 
evaluation of the determinant due to 
computer underflow gr overflow; 
default value is 10 . Indication of 
a need to vary this input is the 
calculation of bending mode shapes 
with frequencies equal to multiples 
of the scan interval (E) 71. 
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El 5 9 Preprocessor Data - (Cont'd) 


(E) Location 

Item 

Description 

1195 

(Control) 

Make nonzero (1.) to activate the 
E159 uncoupled mode preprocessor 
branch of the program. 

1541-1560 

k 

y io 

Thicknesswise mass radii of inertia 
for selected blade breakup, in. Note 
that nonzero values will cause the 
input torsion inertia to be represented 
by both thickness and chordwise radii 
of inertia (input items (R) 121-140 
and (R) 141-160). 

2000 

(Control) 

Case number for output generated by 
E159 branch of the program. 
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Radial Distributions 


(R) Location Item 


Description 


1-20 + Ax 


t 


21-40 

41-60 

t 

M, 

i 

(EIJ f 

61-80 

t 

(ei) e 

O 

o 

1— I 

r — 1 
00 

t 

( I/c) ] 

101-120 

t 

(l/c) ] 

121-140 

t 

k 

y io 

141-160 

t 

k 

z io 

161-180 

t 

S A 

181-200 


y io 

CG 

201-220 


*A 


Nondimens ional blade segment lengths, in 
order from root to tip, maximum of 20 values, 
starting from the offset location. Accuracy 
is generally improved if the last segment 
is small ( < 0.03) . 

2 

Mass of each blade segment, lb-sec /ft. 

2 

Flatwise bending stiffness, lb-in 

2 

Edgewise bending stiffness, lb-in 

Section modulii for flatwise bending, 
root to tip, in^. 

Section modulii for edgewise bending, 
root to tip, in-* . 

Thicknesswise mass radii of gyration of 
blade segments about axis perpendicular 
to chord line and through the reference 
axis, root to tip, nondimens ional with 
respect to R. 

Chordwise mass radii of gyration of blade 
segments about elastic (reference) axis, 
root to tip, nondimensional with respect 
to R. 

Area radii of gyration about elastic axis, 
root to tip, nondimensional with respect 
to R. 

Distances from elastic axis forward to 
airfoil section mass centers, root to tip, 
nondimensional with respect to R. 

Distances from reference (elastic) axis 
forward to edgewise bending neutral axis, 
root to tip, nondimensional with respect 
to R. 
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Radial Distributions (Continued) 


(R) Location 
221-240 


Item 


10 


EA 


241-260 


Ay 


10 


EA 


261-280 


10 


EA 


281-300 


Az 


10 


EA 


301-320 


EB. 


321-340 


EB, 


341-360 


x/M 


Description 

Chordwise distances forward from the 
(extended) root pitch axis to the 
built-in elastic axis (defined herein 
as the locus of section shear centers), 
nondimens ional with respect to R. 

g U xlt-in chordwise elastic axis offset 
change per segment (arc) length distribu- 
tion. This item is a direct statement of 
the sine of the chordwise (forward) 
structural sweep angle distribution; if 
all values of this distribution are input 
as zero, the sweep angle sine distribu- 
tion is computed internally using 
numerical methods from the offset data 
(R)221-(R) 240 , as appropriate. 


Thicknesswise distances from the root 
pitch axis to the built-in elastic 
axis (see above items (R)221 - (R) 240), 
(+) in the normally thrusting direction, 
nondimens ional with respect to R. 

Built-in thicknesswise elastic axis 
offset change per segment (arc) length 
distribution 

Torsional stiffness (to be multiplied by 
twist rate squared), as defined in 
Reference 4, lb-ft^. 


Torsion to edgewise elastic coupling 
stiffness (to be multiplied by twist 
rate), as defined in Reference 4, lb-ft^. 


Constants relating torsional moment to 
torsional stress, root to tip, in”^. 
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Radial Distributions (Continued) 


(R) Location Item 


361-380 



381-400 



401-420 A 0 

B 


421-440 c 


441-460 





Description 

Aerodynamic built-in nonlinear twist 
angle distribution, deg. Since collective 
angle is defined at the 75% span location, 

0 should have a zero value at 75% span. 

B a ... 

Should the structural twist angle distribu- 
tion differ from 0^ , the appropriate data 
must be loaded into a locations (R) 381-(R) 100 ; 

otherwise. 0 will be used for both 
1 B 

aerodynamic a and structural applications. 

Structural built-in nonlinear twist angle 
distribution, root to tip if different from 
aerodynamic twist, deg. See remarks above 
for aerodynamic built-in twist, locations 
(R) 361 - (R) 380 . 

Built-in (structural) twist angle change 
per segment length distribution, root to 
tip deg. Note that this item is a direct 
statement of the built-in twist rate 
distribution, 0^; if all values of this 
distribution are input as zero, the twist 
rate distribution is computed internally 
using numerical methods from the input 
twist angle distributions, locations 
( R) 3 6 1 - (R) 380 or (R) 381- (R) 400 , as 
appropriate . 

Blade chord at center of each segment 
(use for nonconstant chord blades only), 
root to tip, ft. 

Distances from elastic axis forward to 
airfoil quarter chord position, root to 
tip, nondimens ional with respect to R. 

Aerodynamic built-in sweep angle distribution, 
positive aft, deg. 
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Radial Distributions (Cont'd) 


(R) Location Item 


Description 


481-500 



Blade airfoil section skin friction drag 
coefficient , for use with swept airfoil, 
skewed flow option. 


501-520 

v se 


521-540 (Control) 


541-560 (Control) 


Distribution of edgewise nonviscous 
structural damping ( = g / 2 , where g^ 
is the usual structural damping coefficient) 

Distribution of explicit airfoil shape 
similarity index for use with the unsteady 
stalled airloads option, A zero value 
denotes a default quasi-static airloads 
modeling whereas a value of (1., 2., 3.) 
denotes respective airfoil similarity with 
the (NLR-1, NACA 0012 (Mod), Vertol 23010- 
1.58) airfoil sections, all in the mid to 
high subsonic Mach number range. 

Distribution of optional selection modeling 
with the unsteady unstalled airloads . 

A value of (0., 1., 2.) denotes Pad£ 
approximants evaluated using (generalized 
Wagner function, analytic Jordan theory, 
experimental Davis data) aerodynamic sources, 
respectively. 
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Geometric Characteristics 


(G) Location 

Item 

Description 

i 

t 

Rotor tip speed, ft/sec 

2 

t 

R 

Rotor radius, ft. 

0 

b 

Number of blades 

4 

t - 
e 

Nondimensional offset distance of start - 
deformable and/or deflectable portion of 
rotor blade , e/R 

5 

f NSEG 

Number of blade segments used to define 
spanwise variable arrays 

6 

0 

Rotor area solidity (bc/rR) 

7 

0 

1 

equ 

Linear equivalent blade twist angle, 
i.e., difference between tip and root 
built-in angles, positive when tip angle 
is greater (L.E. up) than root angle. 
Note, this input is the default value 
when both the built-in twist angle 
radial distributions have all zero 
values, deg. 

8 

c 

Blade chord if chord is constant, ft. 

10 

B b 

Built-in precone angle, deg. 

11 

<5 

B 

Built-in prelead angle, deg. 
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Dynamic Related Parameters 


Location 

Item 

Description 

i 

NFM 

Number of flatwise bending modes to be 
used (5 max) 

2 

NEM 

Number of edgewise bending modes to be 
used (3 max) 

3 

NTM 

+ 

Number of elastic torsion modes to be 
used (2 max) 

4-8 

i 

to 

w i 

+ 

Flatwise modal frequencies, nondimensional 
with respect to ft, in ascending modal order 

10-12 

l 

(0 

v k 

Edgewise modal frequencies, nondimensional 
with respect to ft, in ascending modal order 

14,15 

+ - 

V 

j 

Torsion modal frequencies, nondimensional 
with respect to ft, in ascending modal order 

30 

(Control) 

In general, a nonzero value invokes the 
hysteretic (nonviscous) structural 
damping formulation in the edgewise 
bending mode equations. Specifically, 
make (1., 2.) to use (the constant 
location (D)31 value, a spanwise 
distribution of values) in this non- 
viscous formulation, respectively. 

31 

se 

Viscous damping equivalent critical 
damping ratio used to approximate 
structural damping in all edgewise 
bending modes 

32 

C w 

se 

Critical (viscous) damping ratio for 
structural damping in flatwise bending mode; 

33 

s 

se 

Critical (viscous) damping ratio for 
structural damping in torsion modes 
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Dynamic Related Parameters (Cont'd 


(D) Location 


Item Description 


34 


39 


K Torsional spring rate connecting root of 

root blade to fixed structure to represent 

control system flexibility, ft-lb/rad. 

A nonzero value will automatically 
introduce the rigid-body feathering 
degr ee-of-f reedom as an additional 
"torsion mode", 

2 

g Acceleration due to gravity, ft/sec ; 

a negative value implies inverted flight. 
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Aerodynamic Characteristics 


(A.) Location 


1 

2 

3 


4 


5 


6 


7 


Item 

P 

a 

oo 

B 


Ac, 


o 


N 

cut-ou 


t 


(c d } 


cut-out 


K 


vim 


Description 

2 4 

Air density, lb-sec /ft 

Speed of sound, ft/sec 

Tip loss, used to define equivalent 
momentum area and three-dimensionality 
corrections to computed two-dimensional 
airloads near the blade tip. 

Increment added to all values of c 

. . d 

obtained from tabulated airfoil data or 
from the analytic NACA 0012 data. 

Airfoil data generally correspond to 
smooth wind tunnel models and A c^ 
is often used to adjust for the 
higher drag of production blades; a 

commonly used value of Acj is 0.002 

o 

Number of blade segments, starting at 
inboard end and defining the cut-out 
region for which the lift and moment 
coefficients are set to zero. 

The drag coefficient used on the first 

N segments, 

cut-out 

Effectivity factor of the induced velocity 
calculated using actuator disk momentum 
considerations in calculating inflow 
angle at a local blade section. Default 
value is 1., corresponding to conventional 
usage of momentum actuator disk inflow. 
This input quantity can be used to 
approximate the effects caused by real 
inflow characteristics as modeled by 
more accurate theories. For such 
usage, the effectivity would typically 
be in the range of 1.0 1,1 


186 


Aerodynamic Characteristics (Cont'd) 


(A) Location 


21 

22 


28 

29 

30-32 


61 

62 

63 


Item 


Description 


V 


Forward flight velocity, kts 


Shaft angle-of-attack measured relative 
to vertical position, deg. Thus, a 
forward thrusting propeller would 
generally have the value -90. 


, 75R 


Blade collective pitch angle as defined at 
the 75% radius, deg. 


X 


Mean rotor inflow ratio 



(Control ) 


(Control ) 


( Control) 


Initial conditions on the "momentum" 
induced velocity components comprising 
a Glauert-like variable inflow description. 
Note that the "vor tic i ty" variable inflow 
(controlled by locations (A)65 and (A)66. 
and the momentum variable inflow can 
be used separately or simultaneously. 

Make nonzero (1.) if airfoil data for a 
nonsymmetric airfoil are to be used. 

Make 1., to invoke the radial flow, swept 
airfoil option. 

Analytic (static) airfoil option. Make 
nonzero (1.) to use the built-in analytic 
approximation to the static NACA 0012. 
airfoil data. 
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Aerodynamic Characteristics (Cont'd) 


Item 

(Control) 


(Control ) 


(Control) 


Description 

Unsteady airfoil data option. A zero 
(default) value invokes the use of the 
conventional quasi-static airloads 
description in both the eigensolution 
and time-history (t-h) solution. 

Nonzero values invoke the following usages 

1. = generalized Wagner function to 

define effective angle-of- 
attack; assumes unstalled aero- 
dynamics, Pade form in eigen- 
solution and tabular airfoil 
data look-up in t-h solu. ; 

2. = UTRC synthesized a, A, a w 

method in the t-h solu, quasi- 
static in eigensolution; assumes 
dynamic stalled aerodynamics 
(see locations (R) 511-540) 

3. = Complete Pade description in 

eigensolution and t-h solu. ; 
t-h solu. uses filtering to 
define constant angles-of- 

attack for use with tabular 

* 

airfoil data look-up; Pade 
airloads are then perturba- 
tional (see locations (R) 541- 
560) . 

4. = Complete Pade description with- 

out filtering in t-h solu. ; non- 
perturbational usage with no 
tabular data look-up. 

Make nonzero (1.) to load (vorticity) 
induced velocity distributions from the 
PANPER code using the input unit code 
specified in input location (S)33. 

Make nonzero (1.) to use the (PANPER) 
induced velocity distributions loaded as 
per input locations (A)65 and (S)33. 



Aerodynamic Characteristics (Cont'd) 


(A) Location 


Item 


68 (Control) 


69 (Control) 


70 (Control) 


71 



Description 

Generalized Glauert (momentum derived) 
variable inflow option. A zero value 
deactivates usage. A value of 1. causes 
the input induced velocity components 
to be used as constants; a value of 2. 
causes the velocity component, v , to 
be varied to satisfy momentum ° 
considerations . 

Input nonzero (1.) to activate use of 
the tabulated time-histories of incremental 
control angle, A 0 


Gust wave option . Make (1., 2.) to use 
the input gust wave function table as 
(a gust factor > (ND) ;an incremental gust 
velocity, fps), respectively. 

Inclination angle of the gust wave, deg, 
A positive value implies an upward 
component . 


72 



Side flow angle of the gust wave, deg. 
A positive value implies a component 
from the starboard side. 


101 


A 0 


.75 


102-149 A 0 75 

251 (Control) 


Number of abscissa-ordinate point pairs 
used to define time-history of A0 (t); 
calculation of this time-history 
is bypassed with a zero value. 

Table of 4 0 ^ abscissa-ordinate pairs; 

U ( A0 ) =' deg.; U (t) - sec. 

Number of abscissa-ordinate point pairs 
used to define the time-history of the 
gust wave function; calculation is bypassed 
with a zero value. 
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Aerodynamic Characteristics (Cont'd) 


(A) Location 


Item 


252-299 f 

gust 


300 


(Control) 


301 


M 


lower 


302 


M 

upper 


305 


(Control) 


306 


(Control) 


307 


308 


309 


(A) 


max 


(a ) 
w 


max 


Description 


Table of gust wave function abscissa-ordinate 
pairs; u(f ) = (see input loc. (A)70); 
U(t) = sec . ® u8t 

Static data utilization option for use with 
the UTRC a-A-a synthesized unsteady 
stalled airloads. A value of (1., 2.) 
invokes the use of (input static data, 
built-in static data) respectively. 

Lower Mach number below which the 
UTRC unsteady airloads method is 
suppressed . 

Upper Mach number above which the UTRC 
unsteady airloads method is suppressed 
(0. - no upper limit), 

Input/Output unit code number for restart 
output of unsteady parameters. ( . 0 = 
input/output is suppressed). 

An input value of (0., 1.) causes the 
UTRC unsteady stalled airloads method to 
(omit, include), respectively, the 
calculation of unsteady drag . 

Maximum (absolute) value of A parameter 
used in UTRC unsteady stalled airloads method 
(0. = no limit) 

Maximum (absolute) value of unsteady decay 
parameter, a , used in unsteady airloads 
formulations (0. - no limit) 

Frequency used for numerically differen- 
tiating inflow angle to calculate A param- 
eter. Typical nonzero values would be 
blade 1st or 2nd torsion mode frequency. 

Zero value invokes numerical differentia- 
tion algorithm based on standard backward 
differencing techniques . Nondimens ional 
with respect to ft. 


190 


Aerodynamic Characteristics (Cont'd) 


(A) Location 
310 


Item 


Description 


(Control) Number of spanwise shape functions 

(Legendre polynomials) used for implemen- 
tation of the Pad£ airloads in the 
eigensolution . (default = 1., maximum 
value = 3 . ) 
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State Variables 


(V) Location Item 


Description 


24 


❖ 


Initial condition on rotor azimuth, deg. 


* 

41-45 q 

w . 
1 


Initial conditions on i 1 th flatwise 
bending mode (nondimensional) rates 


46-48 



Initial conditions on k'th edgewise 
bending mode (nondimensional) rates 


49,50 q Q 

j 


Initial conditions on j ' th torsion 
mode (nondimensional) rates. 


61-65 



Initial conditions on i ' th flatwise 
bending mode deflection. 


66-68 



Initial conditions on k'th edgewise 
bending mode deflections 


69,70 q 

e . 

j 


Initial conditions on j ' th torsion 
mode deflections. 
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Solution Control Parameters 


(S) Location Item 

1 CASE 

6 (Control) 


7 (Control) 

8 A \J/ 


9 


n f 


10 


Description 


Make nonzero (1.) to output print the 
modal integration constants used in the 
eigensolution and, to a limited extent, 
in the time-history solution. 

Make nonzero (1.) to suppress the 
eigensolution 

Azimuth increment used in the numerical 
integration o f the dynamic equations, deg. 
(See section on general information for 
efficient program usage.) 

Number of "flap trials", i.e., maximum 
number of rotor revolutions for which 
the blade time-history will be computed 
in an attempt to obtain convergence to 
periodicity. If a transient response 
is desired for only a portion of one 
rotor revolution the program will 
compute a time-history solution for any 
nonzero fractional value input. 

An identically zero value will cause 
the time-history solution to be 
by-passed entirely. 

Flapping tolerance to within which the 
aeroelastic/dynamic responses must repeat 
on successive revolutions in order for 
the motion to be considered converged to 
periodicity . 
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Solution Control Parameters (Cont'd) 


(S) Location Item 


11 


A ip . 

print 


12 (Control) 


15 (Control) 


21 (Control) 


23 (Control) 


24 (Control) 


25 (Control) 


Description 

Azimuth increment used to present 
printed output of various pertinent 
aerodynamic, dynamic and elastic load 
distributions as well as aeroelastic 
responses and stresses, deg. This 
input quantity should be an integral 
multiple of location(S)8 , deg. 

Make nonzero (1.) if the total (t ransient ) 
time-history is to be output; i.e., 
responses calculated before convergence 
to periodicity is obtained. 

Make greater than zero (1.) for first case 
or when new blade modal data are to be 
input. Program automatically sets this 
control number to (-1.) after each 
loading of modal data. 

Input nonzero (1.) for stress calculations 
using the mode deflection method. Zero 
value defaults to force-integration method 

Make nonzero (1.) if the modal responses 
and hub shears and moments are to be 
(negative) Fourier analyzed after 
periodicity has been obtained. 

Input nonzero (1.) to harmonically analyze 
and output harmonics of flatwise and 
edgewise bending stresses. 

Input nonzero (1.) to harmonically analyze 
and output harmonics of torsional stresses 


26 


Aip 


plot 


Azimuth increment used to form the data 
strings for plotting purposes , deg. 

This input quantity should be an inte- 
gral multiple of location (S) 8. 
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Solution Control Parameters (Cont'd) 


(S) Location Item 


27-31 (Control) 


32 (Control) 


33 (Control) 

34 (Control) 


33 


(Control) 


37 (Control) 


41 


SR 


Description 

Spanwise segment numbers for which 
stress data is to be saved for plotting 
or TSSA purposes (max value = 5.) 

Unit code number of file into which the 
saved time data strings are to be stored 
for subseqeunt plot or TSSA purposes 
(default =12.) 

File unit code number for data transfer 
with PANPER code . 

File unit code number for use in saving end 
conditions for use as initial conditions 
in subsequent (restart) runs and/or with 
the PANPER code. 

Stress calculation supression option . 

Make value equal to (0., 1., 2.) to 
cause suppression of stress calculations 
for (nonoutputted responses only, 
nonoutputted responses and all transient 
responses, all responses), respectively. 

Selection of responses to be used for 
eigensolution linearization. Make 
(0., 1.) to use (initial input modal 
values, modal values azimuthally averaged 
from preceeding run), respectively. 

Sample rate for Transient Spectral 
Stability Analysis (TSSA). (See Reference 32 
for a discussion of this technique.) 

Every (SR)'th point in a transient time- 
history is saved for use in a TSSA. A 
zero value bypasses the TSSA. 


195 



Solution Control Parameters (Cont'd) 


(S) Location Item 


42-44 (Control) 


46 w 

JJ 

47 “u 

48 (Control) 


49 


(Control) 


50 


FREQ 


Description 

Channel selection for each of three 
available for the TSSA. The channels 
available are: 1-5, flatwise bending 
modal responses; 6-8, edgewise bending 
modal responses; 9,10 torsion modal 
responses; 11, 12 & 13, blade tip 
vertical, inplane and torsion deflections, 
respectively. In addition, channels 
14 through 28 are available for stability 
analysis of the blade stresses as selected 
by input locations (S) 27-31. Thus, 
channels 14 thru 18 contain the five 
flatwise stresses, channels 19 thru 23 
contain the five edgewise stresses and 
channels 24 thru 28, the five torsion 
stresses . 

Lower bound of frequency band chosen for 
TSSA nondimensional with respect to ft. 

Upper bound of frequency band chosen for 
TSSA nondimensional with respect to ft. 

Initial estimate of percentage of total 
transient data used in each time dis- 
placed data sample block in TSSA. 

Number of transient ( time displaced ) 
Fourier coefficient calculations made 
to establish modal damping in TSSA; 
maximum value is 200. 

Number of desired resonant frequencies 
to be extracted from frequency band 
defined by input locations (£) 46 and 47. 
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Solution Control Parameters (Cont'd) 



99 (Control) Location used to end a case or series of 

cases. Input (+1.) to end the Loader 
Format data block for the case defined 
by the Loader data and load additional 
cases at the conclusiora of that case. 
Make (-1.) to end the loader data and 
read no further cases. In both 
instances the combined aphameric code 
and word count, ZZ NN (see beginning 
of this section above) should be either 
(Si) or (-1). Note: this entry must 
appear singly on an input card, and 
that card must be the last card for 
the case. 
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III. Blade Mode Shape Data 


included in this data are the radial distributions of the blade (uncoupled) 
flatwise, edgewise and torsion normal mode shapes and their derivatives. These 
quantities are either calculated internally by the (activated) E159 branch of 
the program (and not input as part of the input data) or are explicitly input 
from some previous source. This previous source can be either the output from 
the Ell 59 branch itself (from some previous run), or an equivalent analysis. 

In ejLther case, these data are then input in the following card image format: 


/ NFM 

NEM NTM 

NSEG 

(414) 

cards : 

^F(i) 

F(i+1) 

F(i+2) F(i+3) F ( i+4 ) 

(F18.0, 4F12.0) 


where: NFM, NEM, and NTM are, respectively, the numbers of flatwise bending, 

edgewise bending and torsion normal modes whose mode shapes and derivatives 
are to be input. NSEG is the number of blade spanwise stations for which 
the input modal data are defined. F(i) are the modal functions listed below 
(defined at the i 1 th spanwise stations). Five entries per card are made 
for each F function input for NSEG total entries. The modal functions must 
be loaded in the following order: 
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IV. Components of Variable Inflow 


If location (A) 65 of the Aerodynamics Characteristics Loader Block of 
operational data is nonzero, the following block of variable inflow is 
loaded in: 


card #1: 

^ (VINRAD (I), I = 1,8) 

(8F10.0) 




card //la: 

/ (VINRAD (I) , 1=9, 16) 

(8F10.0) 



card //lb: 

^ (VINRAD (I) , I = 17, . . .) 

(8F10.0) 




card #2: 

^ (VINTAN (I) , 1=1, . . .) 

(8F10.0) 




card //3: 

^ (VINAC (I) , 1=1, . . .) 

(8F10.0) 


VINRAD (I), VINTAN (I), and VINAC (I) are, respectively, the radial, 
tangential, and axial components of the variable inflow at the I f th radial 
station, as provided by the PANPER routine. Card image //la is provided if 
NSEG > 8. Card image //lb is provided if NSEG > 16, etc. Similarly, card 
image //2 is repeated until NSEG tangent values are completed, and card 
image #3 is repeated until NSEG axial values are completed. 


V. Multiple Case Runs 

The above described data setup defines the correct ordering of required 
data blocks for a general G400PROP case, or for the G400PROP portion of a 
more complicated multi-program run stream. When multiple cases are run 
(while remaining within the G400PROP portion of the run stream) the second 
and subsequent cases utilize most of the data input for the first case. 

The following rules apply to the running of mlutiple cases: 
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1. Airfoil data is loaded only for the first case; all subsequent cases 
within the run use the same tabular data, if analytic data is not 
used . 

2. Only those items within th e operational generic (loader) data which 
are to be changed from case to case need to be input. 

3. Item (S) 99 of the operational data controls the running of subse- 
quent cases; a (+1.) value causes a subsequent case to be loaded 
whereas a (-1.) value terminates the computer run after the current 
case. 

4. Unless otherwise specified (by a +1. value for operational data item 
(S) 15) the input modal array data block is used for all cases within 
the run. 

5. Similarly, unless otherwise specified (by a +1, value for operational 
data item (A) 65) and appropriate additional variable inflow data, the 
input variable inflow data block is used for all cases within the run. 

6. Operational data items (S)15 and A(65) discussed above are both auto- 
matically set to zero at the conclusion of the data input for every 
case. 

7. Terminal conditions on the blade azimuth angle, item (V) 24, and on 
the degrees-of-freedom, items (V) 41-50 and (V) 51-70, for any case 

are carried over as initial conditions on these quantities for the sub- 
sequent case. Thus, for some applications, e.g., investigations of 
unstable responses, it would be appropriate to reinitialize these items 
on the subsequent cases. 

When solution flow leaves the G400PR0P portion of a complex run stream, 
the ability to carry over terminal conditions (as initial conditions 
for a subsequent case) and/or any other quantities associated with trim 
is lost. However, a need still exists for preserving these initial 
conditions for subsequent reentries to the G400PROP portions of the run 
stream. As per loader locations (S) 34 (described above) these initial 
conditions are written to and read from the file indicated in this input 
location. 
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PROGRAM OUTPUT DESCRIPTION 


The complete printed output generated by the G400PROP program can be 
classified into the following seven major categories: 

I. Listing of Input Airfoil Data 

XI. Results of Uncoupled Blade Mode Calculation 

III. Listing of Remaining Input Data 

IV. Parameters Calculated from Input Data 

V. Results of Solution Part I - Eigensolutions 

VI. Results of Solution Part II - Time-History Solution 

VII. Results of Solution Part III - Transient Spectral Stability 

Analysis 

This section describes the pertinent output pages associated with each of 
these categories. Generally, output will always be generated only for the 
third category. Output for all other categories depends upon optional 
activation of the uncoupled blade mode calculation, solution parts I, II 
and/or III, and upon optional suppression of the airfoil data. The subsec- 
tions which follow describe, in turn, the details of each of these seven 
categories. Where appropriate, reference is made to descriptions of 
input items described in an above section. 


Listing of Input Airfoil Data 

If static airfoil data is input, then a listing of this data will be 
output for c lt c d and c mc/ , 4 each with the format shown in Sample Page 1. 
First, each of the three aerodynamic section coefficient types is appro- 
priately identified with a label. Next, if multiple spanwise airfoil data 
is input, the radial station at which the data is defined is output. The 
bulk of the remaining data on Sample Page 1 is the actual airfoil data 
wherein each column represents data at one Mach number. Within each 
column, the first line gives the number of angle-of-attack/aerodynamic 
coefficient pairs defining the functionality; the second line is the Mach 
number, and the ensuing line pairs are the angle-of-attack/aerodynamic 
coefficient pairs, where the angles-of-attack are in degrees. This output 
closely follows the input format description given in an above section. 
Note that if the number of radial stations for c^ data sets is input 
negatively, then the total output of the airfoil data is suppressed. 
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Results of Uncoupled Blade Mode Calculation 

A requirement for the use of the basic G400 analysis is a set of 
uncoupled blade modal data consisting primarily of natural frequencies 
and normal mode shapes (and their spanwise derivatives). Since the 
initial development of the G400 code (Reference 1) , the United Technologies 
Corporation program E159 has been added to the G400 code in the form of 
a preprocessor. Since the output from the E159 preprocessor impacts on the 
Loader data used in the G400 proper portion of the code, it must be run 
before the output of the Loader data. Hence, this output category pre- 
ceeds the output of the remaining input data. 

Sample Page 2 and the top of Sample Page 3 echo data input in the "E" 
block of the Generic Loader data. The principal data output on Sample Page 
2 are six pairs of columns defining the distributions of pertinent elasto 
mechanic properties starting at the blade root and progressing outward to 
the blade tip. The lower, major portion of Sample Page 3 consists of four 
columns, the first of which echos the selected blade segment breakup as 
nondimens ionalized by blade radius. The other three columns present at 
each of the selected radial stations, in respective order, the lumped equi- 
valent masses, the effective flatwise and edgewise elastic coefficients 
(spring rates) as calculated from the input stiffnesses and segment lengths. 
These elastic coefficients are used only internally to calculate the natural 
modes. Sample Pages 4 and 5 present, for the selected radial stations, the 
effective values of those elastic stiffness properties which are used by the 
subsequent G400 proper portion of the code. On each of these sample pages, 
the columns give the variations of properties together with the correspond- 
ing selected blade segment breakup. 

Sample Page 6 begins the output of the actual calculations for the 
bending modes. The format of Sample Page 6 is used for both flatwise and 
edgewise modes and the type of modal data being output is clearly indicated. 
The first two rows of Sample Page 6 give the natural frequencies of the 
flatwise modes in increasing order. The units of these frequencies, as 
given in the first row, is either in rad/sec or per rotor rev depending on 
whether the input rotor speed is zero or nonzero, respectively. The second 
row values are always the same respective frequencies, but in units of Hz 
(cycles/sec) . The twelve columns which follow on this sample page consist 
first of descriptors for the selected blade radial stations. The second 
column (labeled X) gives the nondimens ional spanwise stations as measured 
from the axis of rotation, whereas the third column (labeled XFH) gives the 
similar information, but instead from the offset point. The remaining 
three groups of three columns give, for up to three modes, the mode shape 
and derivative information, as indicated. The remaining three columns of 
output at the bottom of the sample page give, for each respective mode, 
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Sample Page 



various internally calculated modal integration constants as indicated by 
the definitions given parenthetically at the left of the page. Note that 
G refers to the respective mode shape array given directly above each of 
these three lower columns. 

Sample Page 7 presents for the selected blade breakup the effective 
values of torsional inertia properties. The second column (labeled IXX) 
presents the section polar torsion inertia, being the sum of chordwise and 
thicknesswise mass moments of inertia. The third column presents the 
difference between the chordwise and thicknesswise mass moments of inertia. 
This difference of inertia distribution is the one which determines the 
so-called "propeller moment" torsional stiffening effect. The remaining 
columns are the mass radius of gyration distributions which, along with 
the output lumped equivalent masses, produce the values in the second 
and third columns. 

Sample Page 8 presents the actual calculations for the torsion modes. 
As with the bending mode output, the first two rows of Sample Page 9 give 
the modal natural frequencies first in either rad/sec or per rotor rev, 
and then in the second row, cycles/sec. The remaining columns follow 
closely the format discussed above for the bending modes, except that only 
the first spanwise derivative is given for each mode. No modal inte- 
grations are performed for the torsion modes. 


Listing of Remaining Input Data 

Output in this category includes a descriptive listing of the Generic 
Loader data and a listing of the components of variable inflow obtained from 
the PANPER code. A description of the Generic Loader data output is 
omitted herein since this output merely duplicates the descriptions already 
given in even greater detail in a previous section. In Sample Page 9 is 
shown the listing of the variable inflow input from a data file (prepared 
by the PANPER code) . This output consists of three columns together with 
a column giving the spanwise station index. The three data columns give 
the cylindrical coordinate system components of the variable inflow. The 
positive directions of the components are as follows: radial (VINRAD) , posi- 

tive outward; tangential (VINTAN) , positive toward the leading edge; axial 
(VINAX) , positive in the normally thrusting position, and all three compon- 
ents have the units of ft/sec. 


Parameters Calculated from the Input Data 

In Sample Page 10 is presented data pertaining to the built-in elastic 
axis offset (structural sweep). The various columns present spanwise 
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INPUTTED VARIABLE INFLOW DISTRIBUTIONS (OUTPUT FROM PROGRAM PANPER) FPS 
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ELASTIC AXIS OFFSET DISTRIBUTIONS 
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Sample Page 10. 


distribution as identified by segment number (N) , and spanwise station (X), 
The third and fourth columns give the edgewise offset information for the 
elastic axis offset and the fifth and sixth colums give the corresponding 
flatwise elastic axis offset information. The quantities Y10EAP and Z10EAP 
are derived quantities and are actually the sines of the structural sweep 
angles in the edgewise and flatwise senses. For small angles, these 
quantities become the slopes (first derivatives) of the offset locations 
(Y10EA and Z10EA) . The quantity DUEAO represents the foreshortening array 
needed to restore the equivalent straight beam back to the originally 
structurally swept position. 

In Sample Page 11 are shown typical modal information for the input 
flatwise and edgewise bending modes. For each such mode, the (nondimen- 
sional) modal frequency and input mode shape and spanwise derivatives are 
listed. In addition, the listing presents the various derived incremental 
deflection correction function vectors which account for blade twists 
(see Eqs. (7) through (14)), and for radial foreshortening. Within a flat- 
wise modal information group, the DVB and DVE arrays correspond to those 
first order Av spanwise functions due to built-in twist and torsional modal 
twist, respectively. The DWWBB, DWWBC and DWWCC arrays correspond to the 
second order AW functions. The various arrays DV2BP, DV2EP, DWW2BBP, 

DW2BCP and DWW2CCP are the first spanwise derivations of the second com- 
ponents (those with superscript "2") of the above discussed arrays, DVB, 

DVE, ..., respectively. The DUEAF deflection array corresponds to the 
bending deflection linear foreshortening accruing from built-in structural 
sweep (see Eq. (30)). The UWE deflection arrays correspond to the flatwise 
bending nonlinear foreshortening (see Eq. (33)). Within the edgewise modal 
information group, the various arrays, DWB, DWE, ... etc., correspond to 
similarly defined spanwise functions involving twist, structural sweep and 
the edgewise modal deflection and spanwise derivative arrays. On Sample 
Page 12 are presented the input torsion modal arrays and natural frequencies 
together with the derived pseudo- torsion mode shape (as defined in Reference 
1) and spanwise derivative. Also presented are the first and second order 
deflection correction functions accruing from structural sweep. 


Results of Solution Part I - Eigensolution 

Sample Pages 13 through 17 present the pertinent details of results 
from the eigensolutions . There are two -distinct eigensolutions calculated 
and the two results are identified, respectively, by the following labels: 

1. LINEARIZED NONLINEAR VACUUM CASE 

2. LINEARIZED NONLINEAR NONVACUUM CASE 
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LINEAR AND NONLINEAR MODAL DEFLECTION VECTORS 


ft xxxxxxx 
a xxxxxxx 

8 aaaaaaa 

g xxxxxxx 

Q XXXXXXX 


xxxxxxx 

XXXXXXX 

aaaaaaa 

XXXXXXX 

'•••••• 

XXXXXXX 


8 aaaaaaa 

> xxxxxxx 

> ••••••• 

a xxxxxxx 


aaaaaaa 

aaaaaaa 

xxxxxxx 

XXXXXXX 


aaaaaaa 

XXXXXXX 

xxxxxxx 


u xxxxxxx 

g aaaaaaa 

> XXXXXXX 
Q ••*•••• 

xxxxxxx 


XXXXXXX 

xxxxxxx 


& xxxxxxx 

aa XXXXXXX 

s gggBsss 

5 xxxxxxx 
xxxxxxx 


******* 

XXXXXXX 


sssssas 

xxxxxxx 

xxxxxxx 


xxxxxxx 
ft xxxxxxx 

g ssssgss 

XX XX XXX 


u xxxxxxx 
m xxxxxxx 

g sssssss 

Q xxxxxxx 
xxxxxxx 


sssssss 

XXXXXXX 

XXXXXXX 


XXXXXXX 

xxxxxxx 

Bssbsss 

xxxxxxx 

XXXXXXX 


XXXxXxx 

xxxxxxx 

XXXXXXX 


XXXXXXX 

Xxxxxxx 

SBsBBSB 

xxxxxxx 

XXXXXXX Q 


S3 ssssgss 

o XXXXXXX 



XXXXXXX 


aaaaaaa 

xxxxxxx 

xxxxxxx 


ft, xxxxxxx 
ft xxxxxxx 

s aaaaaaa 

O XX XX XXX 

xxxxxxx 


pa xXXxxXx ^-n 

g aaaaaaa % 

Q XXXXXXX g 
XXXXXXX Q 


„ aaaaaaa 
g saaaaaa 

XXXXXXX 

• • • • • • • 

XXXXXXX 


xxxxxxx 

aaaaaaa 

xxxxxxx 

xxxxxxx 


xxx X XX x 

ft xxxx xxx 

3 XXXX XXX 

Q XXXXXXX 

xxxxxxx 

xxxxxxx 


xxxx 

H xxx 

aaa 

xxxx 


XXX ft, 


aaaaaaa 

aaaaaaa 

xx xxx xx 
XXXXXXX 


xxxxxxx 

XXXXXXX 

aaaaaaa 

XXXXXXX 

xxxxxxx 


CL XXXXXXX 

ft XXXXXXX 

2 aaaaaaa 

a xxxx xxx 


X X X X > 

xxg^> 
X X X X > 


„ aaaaaaa 
g aaaaaaa 

xxxxxxx 
,*•«••• 
xx xxx XX 


XX^xxxx 

aaaaaaa 

xxxxxxx 

xxxxxxx 


XXXXXXX 

xxxxxxx 


xxxx 

xxxx 

aaaa 

xxxx 


xxx o 
xxx h 

... ft 

xxx ft 


g3 a^saaaa 

g aaaaaaa 

a xxxxxxx 

XXXXXXX 


xxxxxxx 

xxxxxxx 

aaaaaaa 

xxxxxxx 

xxxxxxx 


cl xxxxxxx 
ftl xxxxxxx 

CM XXXXXXX 

s Xxxxxxx 

Q XXXXXXX 

xxxxxxx 


XXXXXXX o 

XXxxxxx 

aaaaaaa g 

XXXXXXX ft 

....... ft 

XXXXXXX ft 


aaaaaaa 

aaaaaaa 

xxxxxxx 

xxxxxxx 


xxxxxxx 

Xxxxxxx 

aaaaaaa 

xxxxxxx 

xxxxxxx 


Xxxx xxx 
xxxxxxx 

XXXX xxx 


xxxxxxx 

xxxxxxx 

aaaaaaa 

xxxxxxx 

xxxxxxx 


aaaaaaa 

aaaaaaa 

x XXXX XX 
X XXXX XX 


xxxxxxx 

xxxxxxx 

aaaaaaa 

xxxxxxx 

xxxxxxx 


ft xxxxxxx 
< xxxXxxx 

g aaaaaaa 

Q xxxx xxx 


xxxxxxx 


xxxxxxx 

xxxxxxx 

aaaaaaa 

xxxxxxx 

xxxxxxx 


aaaaaaa 
x gaaaaaa 

x XXXX XX 


gXgX XXX 

aaaBaBB 

xxxxxxx 

xxxxxxx 


XXXXXXX 

xxxx xxx 

aaaaaaa 

xxxxxxx 


216 



TORSION MODE l MODAL FREQUENCY = .XXXXX 
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In both cases, those nonlinear inertial and elastic terms appearing in 
equations of motions are linearized about the input response variables as 
determined by end conditions from a previous run or by the input initial 
conditions from the "V" generic loader block. In the first case, all aero- 
dynamic terms are omitted; this calculation is for the purpose of obtaining 
the "coupled" modal characteristics. In the second case, the full aero- 
dynamic description is retained; this calculation is for the purpose of 
analysis for aeroelastic stability. In both cases, the eigensolution prob- 
lem is stated in terms of an inertia (A) matrix, a damping (B) matrix, and 
a stiffness (C) matrix, and these matrices are output for both cases. 
Furthermore, for both cases, the extracted eigensolutions (coupled roots 
and frequencies) are listed by root pairs. 

Linearized Nonlinear Vacuum Case 


On the upper half of Sample Page 13 are presented the spanwise distri- 
butions of the major nondimensionalized dynamic and structural properties 
used throughout the program (both for eigensolutions and time— history solu- 
tions). The X and XCEN arrays are the nondimensional distances of the cen- 
ters of the segments from the offset and the rotor axis, respectively. The 
QUAD array constitutes the integration weighting numbers for spanwise inte- 
gration. The THETA- STR array is the pitch angle distribution of the struc- 
tural principal axes and has the units of degrees; in general, it differs 
from the aerodynamic pitch angle distribution appearing on another sample 
page. The two arrays, TWIST-BLT and TWIST-TOT, are the nondimensional 
structural twist rate distributions of the built-in twist and the total 
twist (including initial elastic response), respectively; these arrays have 
the units of radians. The quantities TENSB, EIYB, EIZB and MASSB are, 
respectively, the blade tension, flatwise bending stiffness, edgewise 
bending stiffness and mass distributions, all nondimensionalized by 
appropriate combinations of R, ft and mQ. The (Y10NA)/c and Y10CG) /c 
arrays are, respectively, the distributions of edgewise bending neutral 
axis and mass center offset distances from the elastic axis nondimen- 
sionalized by chord (rather than radius). 

On the bottom half of Sample Page 13 begins the output for the 
linearized nonlinear vacuum case with the three matrices defining the 
eigenproblem. On Sample Page 14 are presented the output formats for the 
cases wherein the roots are either complex conjugates (top half) or a pair 
of real— valued numbers (bottom half). For each root, the eigenvector or 
GENERALIZED MODE SHAPE is presented, normalized to the largest amplitude. 

In the case of a complex conjugate pair, the eigenvector is evaluated using 
only the complex root with the positive imaginary part. The number of 
elements in this vector is identical with the dimension of the A, B, 
and C matrices and represent, in order, the coupled relative responses 
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PART I. E IGENSOLUT IONS OF VARIOUS LINEARIZATIONS OF EQUATION SET - CHARACTERISTIC ROOTS AND COUPLED MODE SHAPES 
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of the (NFM) flatwise modes, (NEM) edgewise modes, and the (NTM) torsion 
modes. The arrays labeled PHYSICAL MODE SHAPE consist of the relative span- 
wise distributions of inplane (Y5) , out-of-plane (Z5) and pitching (THETA) 
deflection components of the coupled mode shape. The quantities labeled 
Y5*, Z5* and THETA* are the velocity or 90 degree out-of-phase components 
of their respective deflection components. The Y5, Y5*, Z5 and Z5* modal 
deflections and velocities are nondimensionalized by R and Q, whereas THETA 
and THETA* have the units of radians. 

Linearized Nonlinear Nonvacuum Case 


Sample Page 15 presents the output format appropriate to the aero- 
dynamically effective eigensolution and is labeled LINEARIZED NONLINEAR 
NONVACUUM CASE. Since there exist two optional forms of static airfoil 
data (analytic NACA 0012 or input tables) as well as two optional forms of 
unsteady aerodynamic methodology available to the eigensolution (quasi- 
static or Pade) an appropriate label is output on this sample page based 
upon the option inputs selected. On the upper half of Sample Page 15 are 
presented the spanwise distributions of the major aerodynamic characteristics 
used to define the pertrubational airloads in the eigenproblem. The X and 
XCEN arrays are defined the same as in Sample Page 13. The units of the 
CHORD array are feet. The angle-of-attack descriptors, THETA-AERO, PHI 
and ALPHA are, respectively, the geometric aerodynamic pitch angle (can be 
different from the structural principal axis pitch angle), the inflow angle, 
and the resulting section angle-of-attack, all in degrees. These angles are 
calculated using the input initial conditions for the response variable 
deflections and velocities. The resulting Mach number (MACH) and section 
coefficients (CL, CD and CM) are used to define the linearization point about 
which perturbational airloads are defined. The quantity KAPPA/U is the span- 
wise variation in aerodynamic moment damping coefficient which, when 
multiplied by the local pitch rate, approximates the potential flow un- 
steady pitching moment coefficient. The quantity (Y10C/4)/C is the span- 
wise distribution of quarter chord offset from the elastic axis nondimen- 
sionalized by chord (not radius). The lower half of Sample Page 15 and all 
of Sample Page 16 present the output format for the A, B and C matrices which 
is similar to the format for the vacuum case matrices. In those cases 
wherein the Pade airloads modeling option is invoked (as indicated in these 
sample pages) , the eigensolution dynamic matrices are augmented to include 
submatrices resulting from perturbational Pade airload augmented state 
variables. These augmented state variables each are governed by their 
own dynamic equations which, together with the differential equations for 
the blade dynamics, define submatrix partitioning within each of the now 
augmented A, B and C matrices. Thus, the two A matrices shown in Sample 
Page 15 are, respectively, the upper left and lower left submatrices com- 
prising the augmented A matrix. Sample Page 16 contains the similar outputs 
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Sampl e Page 17. 


for the augmented B and C matrices. The two B matrices are, respectively, 
the upper left and lower left submatrices comprising the augmented B matrix. 
The two C matrices are, respectively, the upper left and upper right sub- 
matrices comprising the augmented C matrix. With the exception of the lower 
right C submatrix, which is always a diagonal matrix, all other submatrices 
are null and their output is suppressed. The eigenvalues calculated for 
the nonvacuum cases are output in the same format as is described above for 
Sample Page 14. 

Should one of the roots in the nonvacuum eigensolution be unstable, as 
indicated by a positive root or real part of a complex pair, an output 
listing of the force phasing matrices appropriate to the instability is 
generated and outputted as depicted in Sample Page 17. These matrices, 
having the same size as the A, B, and C dynamic matrices, enable the various 
destabilizing forces to be identified; descriptive material for their 
definition and interpretation are contained in Reference 33. 


Results of Solution Part II - Time-History Solutions 

Sample Pages 18 through 22 present the pertinent details of the results 
from the time-history solutions (i.e., transient aeroelastic responses). 

The first row of parameters following the page title represents, for the sub- 
sequent time-history solution, the parameters defining the flight condition. 
These consist of the various control angles (in degrees) , the inflow and ad- 
vance ratios, LAMBDA and MU, respectively, and the initial nondimensionalized 
values of the "momentum 11 induced velocity components. The remainder of 
Sample Page 18 comprises the typical azimuthal listing; this listing is output 
for every azimuth angle which is a multiple of the print azimuth increment, 
loader input (S)ll. 

Azimuthal Printouts 


The first line appearing on all subsequent azimuthal printouts gives 
the rotor azimuth angle, revolution number and time. In addition, the three 
components of the instantaneous Glauert (momentum induced) inflow are out- 
put. The remainder of the azimuthal printout consists of 3 main groups of 
result quantities. The first of the three groups on this sample page lists 
the spanwise distributions of the pertinent aerodynamic quantities, the 
format of which depends on the choice of unsteady airloads modeling selected 
(input location (S)64). 

Uns^t^ady^Mrl^oad^ Ut^l^zin^ an_E£u_iva.l_enjt An^l^-o^f^A^tta-ck 

Unsteady airloads options wherein the airloads depend on a single 
"effective" angle-of-attack are: i) the quasi-static modeling ((A) 64 = 0.), 
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PART II. TIME HISTORY SOLUTION OF COMPLETE (NONLINEAR) EQUATION SET - AEROELASTIC TRANSIENT RESPONSES 
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ii) the unstalled generalized Wagner functions modeling ((A)64 = 1.), and 

iii) the UTRC synthesized unsteady stalled modeling ((A)64 = 2 .)* For 

these cases, the first seven aerodynamic quantities are: LAMBDAS, GAMMA, 

PHI, ALPHA, MACH, A, and ALPHAW. These quantities are, respectively, the 
aerodynamic sweep angle due only to radial flow (in degrees), the total 
aerodynamic sweep angle including structural sweep (in degrees), the in- 
flow angle and angle-of-attack (both in degrees), Mach number, nondimen- 
sional angle-of-attack rate, and the unsteady decay parameter (in degrees). 

Uns^t_eady_Ai L rl_ o a.ds^ UtjLljLzijij* the ^^e_F£rmuj L a_tio^n_ 

i 

The Pade airload options inherently assume that the pitch and plunge 
of the airfoil are independent variables defining the airloads. For either 
of the Pade options ((A)64 = 3. or 4.), the first seven aerodynamic quan- 
tities are those shown in Sample Page 18: PHTILDA, THTIL, PHI, ALPHA, MACH, 

DCL and DCM. The first four of these quantities are, respectively, the per- 
turbational inflow and pitch angles away from the mean values defining 
the mean angle-of-attack, the total inflow angle, and the total angle-of- 
attack, all in degrees. The remaining three items are, respectively, the 
Mach number, and the per trubational unsteady lift and moment coefficients. 

If input item (A)64 is equal to 3., the analysis filters out the steady 
(low frequency) parts of the inflow and pitch angles so as to use the Pade 
airloads on only the perturbational airloads (static airloads on the filtered 
low frequency angles). If input item (A) 64 is equal to 4 . , the Pade air- 
loads description is used on the total inflow and pitch angles and the per- 
turbational coefficients become the total coefficients. 

Remainder of Azimuthal Printout 


The aerodynamic coefficients CL, CD and CM are self-explanatory and non- 
dimensional. The airload distributions in the Z 5 and y^ directions, SAZ5 
and SAY5 , respectively, have the units of lb/in., the aerodynamic pitching 
moments about the x^ and y^ axes, MAX5 and MAY 5 , respectively, have the units 
of lb-in. /in. The second of the three groups on Sample Page 18 lists per- 
tinent spanwise distributions of a structural dynamic nature. The vertical 
and inplane deflections are those, respectively, in the Z 5 and y,- directions 
and have the units of in. The torsional deflection has the units of deg. 

The quantities SDZ5, SDY3 and MDX5 are "semi- dynamic” load distributions. 

These distributions are dimensionally the same as those resulting from aero- 
dynamics, but arise instead from all the dynamic effects except the doubly 
time differentiated ones. The quantity MEX9 is the nonlinear elastic torsion 
moment distribution as calculated using the AEI implementation (Eq. (38)); 
it too has the units of lb-in. /in. All stress quantities have the units of 
lb/in. 2, whereas the torsion moment has the units of lb-in. The third of 
the three groups on Sample Page 18 lists the modal responses, their non- 
dimensional time derivatives and "right-hand-side" excitations. Specifically, 
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for each modal response variable (column) are given the instantaneous 
displacement, velocity, acceleration and generalized excitation (elements 
on right-hand side of the modal equation) . 

After the time-history solution has either converged to periodicity 
or run to maximum flapping trials (input location (S)9), various integrated 
loads are calculated for one final blade revolution to form the aero- 
dynamic performance and stress results depicted in Sample Page 19. For 
each of eight performance quantities, results are presented in nondimen- 
sional coefficient form, in nondimensional form divided by solidity, and 
in actual dimensional form. Note that ten dimensional quantities are listed 
and the units are lb for forces and lb-ft for moments, as appropriate. 

The quantity EQU. DRAG (lb) represents the combined power expended by the 
rotor due to rotor rotation (torque) and translation (drag) divided by 
flight speed. 

The next line duplicates the parameters defining the flight condition 
and includes additional quantities which depend on the integrated performance 
for evaluation. At the beginning of the time-history calculation, it is not 
known which part of the inflow ratio being used is due to ram effects and 
which due to momentum induced effects. Once the integrated rotor thrust 
is calculated, however, the induced portion of the inflow can then be cal- 
culated using the simple usual momentum formula derived for flight in an 
infinite continuum. The complementary portion of the inflow represents 
the ram effect from which the shaft angle-of-attack ALPHA S, in degrees, 
can be calculated. The quantity VEL ACT is the actual forward flight 
velocity, in knots, consistent with the advance ratio used and the shaft 
angle-of-attack. For finite forward flight speeds EQU. L/D is the lift 
divided by the equivalent drag; for hovering cases this quantity is the 
figure of merit. PAR. AREA, the rotor parasite (drag) area, in square feet 
is the rotor drag divided by dynamic pressure. The control angles, A1S, BIS, 
THETA 75 and the shaft angle-of-attack all have the units of degrees. The 
power absorbed by the propeller from the airstream in kilowatts is given 
by the quantity KWATT. It will always be of opposite sign from the horse- 
power. The remainder of Sample Page 19 consists of reductions of the 
various stresses given in the azimuthal printout in terms of median and 
1/2 peak-to-peak values. 

Once the time-history solution has converged to periodicity, the pro- 
gram optionally performs harmonic analyses of the azimuthal variations of 
various response quantitites. The outputs of these harmonic analyses are 
depicted in Sample Pages 20 through 22. In each of these sample pages, the 
harmonic information for each response variable is contained in the appro- 
priate horizontal band of five rows. The harmonics are listed by columns 
up to a maximum of 10 harmonics. All harmonic analysis output depicted on 
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HARMONIC ANALYSIS OF BLADE RESPONSES 
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Sample Page 20. 
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Sampl e Page 21. 
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Sample Page 22. 


these sample pages assume a negative harmonic content form in keeping with 
the (negative) harmonic form conventionally assumed for the blade pitch 
control and rigid flapping angles. For each harmonic of response variable 
five quantities are output; these quantities are, respectively, the (nega- 
tive) cosine and sine components, the equivalent amplitude and phase angle, 
and lastly, the amplitude of the harmonic relative to all the other har- 
monic amplitudes output. Sample Page 20 depicts the harmonic analyses of tie 
dimensionless modal response variables selected wherein QW(I) , QV(K) and 
QT(J) are, respectively, the (I) flatwise, (K) edgewise and (J) torsional 
uncoupled mode responses. 

Sample Page 21 depicts the harmonic analyses of the total shears and 
moments exerted by one blade to the hub. In contrast to the steady hub 
loads listed in the AERODYNAMIC PERFORMANCE AND STRESSES output (Sample 
Page 19) which are calculated by integrating only the aerodynamic load dis- 
tributions, the total hub loads, which are herein harmonically analyzed, are 
calculated by similarly integrating the combined aerodynamic and the dynamic 
load distributions. The longitudinal, lateral and vertical hub shears com- 
prising the first three quantities of this sample page all have the dimen- 
sions of lb and are defined in the x^-(aft), y^- (starboard) , and z^- (up and 
along axis of rotation) axis directions, respectively. The roll, pitch 
and yaw moments comprising the latter three quantities on this sample page 
have the dimensions of lb-ft and are defined positive (using the right-hand 
rule) about the x^-, y^-, and z^-axes, respectively. Note that the aero- 
dynamic rolling moment whose output is depicted in Sample Page 19 is defined 
positive starboard side down and is opposite from the harmonically analyzed 
total rolling moment depicted in Sample Page 21. Sample Page 22 depicts the 
harmonic analysis of the flatwise stresses at the center of each of the span- 
wise segments. A similar output listing is provided for both edgewise and 
torsional stresses. 


Results of Solution Part III - Transient 
Spectral Stability Analysis 

Transient time-history solutions are often difficult to interpret for 
quantitative stability information. This is due to the fact that the total 
responses so calculated inherently consist of several component modes simul- 
taneously and transiently approaching (or departing from) multi-harmonic 
periodicity with a wide range of natural frequencies and inherent damping 
levels. The extraction of the component responses at discrete frequencies 
in order to examine their individual attenuation characteristics is the pur- 
pose of the Transient Spectral Stability Analysis (TSSA) portion of Program 
G400PROP. The details of this analysis, which utilizes Fourier Transform 
techniques, are beyond the scope of this report but are treated in Reference 
32. 
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Essentially the TSSA first performs Fourier transformations of selected 
time-history data strings, which have been previously generated in the time- 
history solution portion of the analysis (Solution Part II) and appropriately 
saved. The purpose of the Fourier Transform is to identify, within these 
time-histories, those frequencies whose amplitudes are relatively largest 
and which are herein denoted as "resonances”. The TSSA then calculates the 
transient behavior of the extracted amplitudes of these resonances over the 
time-history time interval and estimates equivalent linear stability indices 
(characteristic exponent, critical damping ratio, and time to half amplitude) 

Sample Pages 23 and 24 depict the output typically generated by the 
TSSA. The sequence of output depicted is duplicated for each of the tran- 
sient response channels selected. Sample Page 23 depicts the output 
generated by the Fourier Transform frequency indentif ication portion of the 
TSSA. Shown at the top of the page is the transient response channel being 
analyzed and the frequency range wherein resonance identification is desired 
(input locations (S)46 and 47). The series of five output items to follow 
consist of parameters defining the numerical Fourier transform; note that 
the results of the TSSA incorporate a time nondimensionalization based on 
rotor speed, ft. The tabulation of the Fourier Transform follows wherein, 
for each frequency (harmonic of the fundamental as determined by the total 
nondimens ional time interval), the real and imaginary parts, the square of 
the amplitude and the logarithm to the base 10 of the amplitude are out- 
putted. Generally, this tabulation will consume more than the one page 
indicated in Sample Page 23. After this listing is completed, those fre- 
quencies and their respective square amplitudes which are found to be 
resonances, as defined above, are listed. 

Aside from their several nonlinearities, the dynamic equations of 
motion for propeller blades implicitly contain several linear aerodynamic 
terms which, under conditions of yawed flight, can become periodic. It is 
not expected then, that the aeroelastic time-history responses generated 
by these equations should manifest Floquet Theory characteristics (see 
Reference 34). In particular, the Fourier Transform is capable of iden- 
tifying "multiple resonances" which are separated by (plus or minus) 
multiples of the rotor frequency and which would be found to have approx- 
imately the same damping level as measured by characteristic exponent. 

Hence, the resonant frequencies found by the resonance search are further 
screened to extract only those frequencies with distinct noninteger values 
and which, within the set having the same noninteger values, have the largest 
transform magnitudes. These extracted frequencies are herein denoted "funda- 
mental resonances" and are the only ones examined further for stability in 
the TSSA. 
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Sample Page 24 depicts the results of frequency fine tuning and 
response stability estimation for each of the fundamental resonances ex- 
tracted earlier in the TSSA. The results for each of these frequencies are 
presented in columnar fashion. The top horizontal blocks of output represent 
the frequency fine-tuning results. Of most practical importance are the 
values labeled OPTIMIZED FREQUENCY which are, in nondimens ional (per rotor 
rev) form, the best estimates of the frequency of the fundamental resonant 
frequencies. These frequencies are obtained by an optimization technique, 
the details of which are beyond the scope of this report. The remainder 
of the output depicted on this sample page (for each fundamental resonance) 
consists of three horizontal blocks of output representing various estimates 
of the effective damping characteristics. These three types of blocks are 
best explained by considering, for each of the fundamental responses indicated 
in Sample Page 24, the variation in the natural logarithm of the magnitude 
of resonant frequency content with (nondimensional) time. If these ampli- 
tude logarithms attenuate with time, then the response with that frequency 
content (mode) is deemed stable, and conversely the slope of that attenua- 
tion with time is a measure of the effective linear damping; in the analysis 
this slope is obtained by a simple least-square fit. It may happen that the 
variation of amplitude logarithms with time is neither monotonic increasing 
nor decreasing in which case a condition of maximum or minimum amplitude 
is defined. By weighting the least-square fit either uniformly or with an 
appropriate function accentuating the initial or terminal ends of the ampli- 
tude logarithms data string, the three latter horizontal blocks of output 
depicted in Sample Page 24 are generated. Within each of these blocks, the 
first quantity depicted is the nondimensional CHARACTERISTIC EXPONENT, which 
is analogous to and interpreted in the same way as the real part of the 
eigenvalue discussed in the output for Solution Part I. The REVS to (MAX/ 

MIN) AMPL is an indication of the asymptotic behavior of the component re- 
sponse. STANDARD DEVIATION is the root-mean-squared error achieved in the 
least-square curve-fit and is an indication of the regularity of the ampli- 
tude logarithm function and of the accuracy of the stability estimation. 

Based upon the OPTIMIZED FREQUENCY outputted at the top of the sample page, 
the equivalent CRITICAL DAMPING RATIO is calculated from the characteristic 
exponent using standard formulas. Finally, the output item labeled REVS TO 
HALF AMPLITUDE is the third alternate way in which the equivalent linear 
damping result is presented. 

In Sample Page 25 is depicted the typical additional page of output 
generated at the beginning of every case following the first case of a 
multiple case sum. The two columns depict, respectively, the location num- 
bers and data values for the newly input data distinguishing the present 
case from the previous one. This feature is intended solely as an ease of 
usage output to assist in data management. 
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Sample Page 24. 



PROGRAM G400PROP- - CASE 2 OF A MULTI-CASE RUN. INPUT DATA REPRESENTING CHANGES TO THE PREVIOUS CASE ARE AS FOLLOWS: 


■■"aHaHHHHHHHHHHg 

^xxxxxxxxxxxxx 

XXXXXXXXXXXXXXXX 


^XXXXXXXXX 



xxxxxx 
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PROGRAM BACKGROUND MATERIAL 


General Information to Facilitate Operation 
of Program and Improve Efficiency 


Aside from considerations of the details of the aeroelastic modeling 
which are covered in various above sections and in Reference 1 and 3, ad- 
ditional attention should be paid to maximizing both the efficiency and 
accuracy of the implemented numerical solutions of the dynamic equations* 

This subsection presents material describing areas of concern with respect 
to these numerical methods, and more importantly, ways of dealing with them 
by proper input procedures. 

Blade Segment Selection 

Two decisions must be made in selecting a proper distribution of blade 
segment lengths: how many segments should be used, and where segments 

should be either sparsely or densely packed. The G400PROP code incorporates 
a maximum number of twenty segments, up from the maximum of fifteen offered 
in the earlier versions of G400 (References 1 and 3) . Capability to use 
twenty segments should not be confused with a general need to use all this 
capability in every application. 

Various criteria can be used to guide the program user in making an 
efficient blade breakup selection: 

!• Generally any one segment should not exceed 15 percent of the span. 
This criterion is subjective in that it is based on accumulated 
user experience. 

2. The segment density should be greatest at the innermost portion of 
the blade for the E159 part of the program (uncoupled mode calcula- 
tions) and at outermost portion for the G400 proper part of the pro- 
gram. The requirement for greater blade detail at the root blade 
portion in the E159 calculation stems from the fact that here the 
elastic strain energy is most heavily concentrated and has the most 
variability. It follows that accurate modeling of the equivalent 
springs used in E159 is enhanced by a finer breakup here. The re- 
quirement for greater blade detail in the blade tip portion in the 
G400 proper calculations stems from the concentration here of in- 
ertial and aerodynamic loadings. The aerodynamic loads are espe- 
cially subject to greatest variability at the tip sections. 
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The segment density should be also guided by the specific details 
of the blade in question. Any blade portion which has locally 
concentrated properties should have greater segment density. Also, 
segment boundaries should be selected to conform to the geometry 
inherent in the blade planform. 

4. With some initial extra effort in the preparatory stage, an effi- 
cient breakup can be used wherein the need for a dense breakup at 
the root can be relaxed in the G400 proper portion of the code. 

This effort consists of first running the E159 preprocessor sep- 
arately with a dense breakup at the root to maximize the accuracies 
of the natural frequencies. Then, various inboard segements are 
selectively eliminated or modified from the mode shapes and other 
distributed data before subsequent input to G400 proper part of 
the code. In this manner the elastic modeling accuracy is preserved 
(through retention of the accurate natural frequencies) while re- 
ducing the all over segment count used in the more expensive G400 
proper aeroelastic calculation. 

Input of Differentiated Data 

The coordinate transformations formulated for G400 require two sets of 
spanwise differentiated data which must be explicitly gleaned from the 
geometry of the blade design: structural twist rate and structural sweep 

rate. Although the G400 codes provide for internal numerical differen- 
tiations of these quantities, actual designs often include abrupt spanwise 
variation which cannot be so differentiated efficiently. Consequently, the 
G400 input list includes a direct input of rate related data, and use of 
this input is generally recommended for increased accuracy. 


The method selected for input of rate data on these items is based on 
the assumption that the rates are constant over their respective segment 
lengths. To make the input numbers more meaningful and to minimize data 
preparation calculations by the user, the rate data are input as respective 
changes in the variables (either twist angle and/or elastic axis offset) over 
each segment length. Thus, the actual derivatives are calculated internally 
by division by each segment length and the user is freed of this chore. One 
advantage of this input format is that the resulting numerical values input 
provides a quick check of the data. All such changes can be easily summed 
to yield the integrated change over the whole blade, for comparison with 
input root to tip values of the variables themselves. 
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Temporal Numerical Integration 


As discussed in Reference 1, temporal integration of the higher differ- 
entiated response variables to obtain the lower ones is achieved in the 
G400PROP program using a variant of the Adams integration algorithm. The 
selected algorithm is defined by means of the azimuthal integration step 
size, Vip , and the integration frequency, oi. 

The integration step size should be an integral divisor of 360; a proper 
choice depends on the maximum coupled frequency inherent in the various aero- 
elastic responses. A reasonable upper limit for Vip is 30 divided by the max- 
imum such frequency in per rev. Values of Vip greater than this upper limit 
will compromise the integration accuracy and, for sufficiently large values, 
will cause the computed responses to develop "numerical" instabilities. As 
a corollary, a check on any response which is predicted to be unstable by 
the analysis, is to rerun the case with a reduced integration step size to 
test for the possibility of the unstable response being merely a numerical 
instability. 

For each response degree-of-freedom a different integration frequency, 

“» is used in the integration algorithm; this frequency is, for each of the 
elastic modes, the respective input natural frequencies (locations (D)4-8, 
(D)10-12, and (D)14,15. In addition to defining modal stiffnesses and in- 
tegration frequencies, the input frequencies serve yet another purpose. As 
noted above, the proper value of integration step size, V\p, varies inversely 
with the maximum modal frequency. Thus, run times (caused by reduced step 
size) will significantly increase as any one modal frequency increases. 

Since any degree-of-freedom exhibiting a large natural frequency tends to 
respond quasi-statically, i.e., as if the acceleration (*q) term were neg- 
ligible, a reasonable approximation to the response calculation is to avoid 
the numerical integration of the *q[ term entirely and treat the response 
quasi-statically. This option can be invoked for any such high frequency 
mode by input of a negative frequency; a negative sign will not affect the 
proper usage of the frequency in the calculation of the dynamic equations. 

Note that this optional response calculation can be invoked singly or in 
combination for any of the elastic modal responses (negative values in any 
of locations (D)4-8, 10-12, and 14-15). 

State Vector Initial Conditions 

For many applications useful results can be obtained from the G400PROP 
code with little or no attention paid to the input of initial conditions 
(all elements of the (V) Loader block and locations 30-32 of the (A) Loader 
block). Three situations exist, however, wherein appropriate and accurate 
initial conditions should be input to maximize the usefulness of the analysis. 
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Stability £a^c»il£t^ons_f^r_C£n^i^u£a_ti£ns With_S_tr£n& Nonlinea£i_ti£s_ 

The literature abounds with examples of indentif iable classes of 
rotor instabilities which involve some form of essential nonlinearity. 

The inclusion of such nonlinear characteristics in a comprehensive aero- 
elastic analysis such as G400PROP typically requires making the stability 
calculations for the blade in its predeflected equilibrium position. This 
is true for some forms of instability which are best analyzed using the 
time-history solution, and it is especially true for those for which the 
eigensolution is most appropriate. Because the time-history solution 
inherently includes modelings of all excitations it is the logical tool to 
use for obtaining these equilibrium conditions. The principal difficulty 
with this approach is that these equilibrium conditions are governed by 
the same equations of motion which govern the transients, including any 
instability should it exist. Hence, the presence of an inherent instability 
can interfere with any attempt to obtain the equilibrium conditions, and some 
form of filtering, to suppress the unstable transient, is therefore required. 
The G400PROP code does not have an explicitly dedicated filtering algorit m 
for this purpose, however, and some other ad hoc measure must be used. Two 
such measures which could be used are briefly described below: 

1. Use of exaggerated amounts of structural damping ((D) loader block 
locations 31-33) in a precalculation run using the time-history 
solution. By artificially increasing these damping levels by at 
least one order of magnitude the inherent instability transient 
responses can be quenched while leaving the final equilibrium point 
solution unaffected. The end conditions from this precalculation 
run define the equilibrium initial conditions for a subsequent 
stability run which would, of course, include the restored values 
of structural damping. 

2. Use of the quasi-static solution on one of the degrees-of-freedom 
critical to the instability in a precalculation run (see the above 
subsection) . The time-history solution would then be used in a 
manner similar to the above described method. 

Cal : ibrat£d_Exci L tatJ.ons_of^ Tran£ients_ 

For those cases wherein the instabilities are best investigated using 
the time-history solution, the initial conditions provide a convenient method 
for exciting the transient responses in an unambiguous and calibrated manner. 
This can be accomplished by selecting a critical degree-of-freedom and as- 
signing to it a rate initial condition ((V) loader block locations 41-50) 
equal to its (nondimens ional) natural frequency times an appropriate amp 
litude (in radians) . 
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Restart Calculations 


For some cases involving time-history calculations insufficient rotor 
revolutions may have been selected to define the dynamic phenomenon under 
study and a continuation of the case is required. To this end use should be 
made of the built-in feature of the G400PR0P code to record end conditions 
on the control angles, components of induced velocity and the blade deflec- 
tion state vector. The code provides for automatic output of punched card 
images at the end of the run and, if Loader location (S) 34 is input with 
a nonzero value appropriate to an available file unit number, to that file 
as well. Note, however, that these end conditions will be so recorded only 
if the run makes a normal completion. Premature aborting of the run because 
of excessive response amplitudes, for example, will suppress this output. 
Finally, once the end conditions are recorded they can then be used as 
initial conditions for subsequent runs. 

Computer Run Times and Storage Requirements 
Typical Run Times 


The run times required for making typical calculations with the G400PROP 
code depend on two main factors: the computer type and the extent of the 

aeroelastic problem being solved. This subsection provides information which 
should be of help in bracketing the run time to be expected in any specific 
application. Clearly the run times are in proportion to the degree of so- 
phistication selected from what is avilable in the code. For example, use 
of analytic static airfoil data, quasi-static airloads and reduced numbers 
of both modes and radial stations will significantly reduce the run time, 
albeit at an expense in accuracy. The data presented in Table X was obtained 
from running various problem configurations with the code on the UTRC UNIVAC 
1110/81A. Usage of the program with complete variation of all possible com- 
binations of options is clearly impractical since it would require a prohib- 
itively large number of runs. The cases selected for inclusion in Table X, 
however, should help in identifying the costs of major calculation features. 
The following points can be made from the results of Table X: 


^ * Eigensolutions are roughly an order of magnitude cheaper than time- 
history solutions. 

2. The inclusion of Pade airloads increases the run time of eigensolu- 
tions roughly 50%. 

3* Significant reductions of time are obtainable by proper tailoring 
of the run to only those capabilities needed. 
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TABLE 



Fade approximant 
unsteady stalled 




4 . 


The additional cost of running the transient spectral stability 
analysis (TSSA) is minimal. For the TSSA activated sample case 
10, 256 point time series were used with a moving block size of 
60% and the use of' 100 block displacement calculations. 

Storage Requirements 

The storage requirements for the G400PROP code as installed on the 
UNIVAC 1110/81A are as follows: total IBANK (instruction required storage) 

is approximately 75K (decimal) locations and the DBANK (data required storage) 
is approximately 146K locations for a total of ^221K. The program lends it- 
self reasonably well to overlay techniques and these storage requirements 
have been reduced to ^45K IBANK and 105K DBANK for a total of ^150K on the 
UNIVAC. A major identifiable driver on these storage requirements is the 
eigensolution capability which requires 46K additional directly identifiable 
storage in the overlay structure. For purposes of comparison with other 
computers, G400PROP was written and compiled in FORTRAN 77, and one location 
on the UNIVAC 1110/81A consists of a 36 bit word. 



REFERENCES 


1. Bielawa, R. L. : Aeroelastic for Helicopter Rotor Blades with Time- 
Variable, Nonlinear Structural Twist and Multiple Structural Redundancy- 
Mathematical Derivation and Program User's Manual. NASA CR-2638, 

February 1976. 

2. Bielawa, R. L. : A Second Order Nonlinear Theory of the Aeroelastic 
Properties of Helicopter Rotor Blades in Forward Flight. Ph.D. Thesis, 
Massachusetts Institute of Technology, June 1965. 

3. Bielawa, R. L. : Aeroelastic Analysis for Helicopter Rotors With Blade 
Appended Pendulum Vibration Absorbers - Mathematical Derivations and 
Program User's Manual. NASA CR-165896, June 1982. 

4. Houbolt, J. C., and G. W. Brooks: Differential Equations of Motion for 
Combined Flapwise Bending, Chordwise Bending and Torsion of Twisted Non- 
uniform Rotor Blades. NACA Report 1346, October 1956. 

5. Hodges, D. H. and E. H. Dowell: Nonlinear Equations of Motion for the 
Elastic Bending and Torsion of Twisted Nonuniform Rotor Blades. NASA 
TN D-7818, December 1974. 

6. Rosen, A. and P. P. Friedmann: Nonlinear Equations of Equilibrium for 
Elastic Helicopter or Wind Turbine Blades Undergoing Moderate Deformation. 
UCLA School of Engineering and Applied Science Dpt. UCLA-ENG-7718 (NASA 
CR-159478), December 1978. 

7. Rosen, A.: The Effect of Initial Twist on the Torsional Rigidity of 
Beams - Another Point of View, ASME Journal of Applied Mechanics, Vol. 

47, 1980, pp. 389-392. 

8. Hodges, D. H. : Torsion of Pretwisted Beams due to Axial Loading, ASME 
Journal of Applied Mechanics, Vol. 47, 1980, pp. 393-397. 

9. Shield, R. T. : Extension and Torsion of Elastic Bars With Initial Twist, 
ASME Journal of Applied Mechanics, Vol. 49, 1982, pp. 779-786. 

10. Bisplinghoff , R. H. , H. Ashley, and R. L. Halfman: Aeroelasticity . 

Addison - Wesley Publishing Company, Inc., 1957. 

11. Bielawa, R. L.: Blade Stress Calculations - Mode Deflection vs. Force 
Integration, Journal of the American Helicopter Society, Vol. 24, No. 3, 
July 1978. 


246 



12. Gangwani, S. T. : Prediction of Dynamic Stall and Unsteady Airloads on 
Rotor Blades. Journal of the American Helicopter Society, Vol. 27, No. 4, 
October 1982, pp, 57-64. Also Preprint No. 81-1, AHS, 37th Annual 
National Forum, New Orleans, LA., May 1981. 

13. Dadone, L. U.: Two-Dimensional Wind Tunnel Test of an Oscillating Rotor 
Airfoil, Vol. II. NASA CR-2915, December 1977. 

14. Gray, L. and J. Liiva: Two-Dimensional Tests of Airfoils Oscillating 
Near Stall, Vol. II, Data Report. USAAVLAB TR-68-13B, USAAMRDL, 

Ft. Eustis, VA., April 1968. 

15. Beddoes, T. S.: A Synthesis of Unsteady Aerodynamic Effects Including 
Stall Hysteresis. Vertica, Vol. 1, pp. 113-123, 1976. 

16. Carlson, R. G., G. L. Commerf ord , and P. H. Mirick: Dynamic Stall Mod- 
eling and Correlation with Experimental Data on Airfoils and Rotors. 
Presented at the AHS/NASA Specialists* Meeting on Rotorcraft Dynamics, 

NASA SP-352, Rotorcraft Dynamics, Paper 2, February 1979. 

17. Gangwani, S. T. : Synthesized Airfoil Data Method for Prediction of 
Dynamic Stall and Unsteady Airloads. NASA CR-3672, February 1983. 

18. Jordan, P. F.: Aerodynamic Flutter Coefficients for Subsonic, Sonic and 
Supersonic Flow (linear Two Dimensional Theory), RAE RM No. 2932, April 
1953. 

19. Vepa, R.: Finite State Modeling of Aeroelastic Systems, NASA CR-2779, 
February 1977. 

20. Ballhaus, W. F. And P. M. Goorjian: Implicit Finite Difference Computa- 
tions of Unsteady Transonic Flows About Airfoils, AIAA Journal, Vol. 15, 

No. 12, December 1977, pp. 1728-1735. 

21. Ballhaus, W. F. and P. M. Goorjian: Efficient Solution of Unsteady Tran- 
sonic Flows About Airfoils, Paper No. 14, AGARD Conference Proceedings 
No. 226, Unsteady Airload in Separated and Transonic Flows, 1978. 

22. Hildebrand, F. B.: Introduction to Numerical Analysis. McGraw-Hill Book 
Co., Inc., New York, N. Y., 1956. 

23. Davis, S. D. and G. N. Malcolm: Experimental Unsteady Aerodynamics of 
Conventional and Supercritical Airfoils, NASA TM J1221, August 1980. 


247 



24. Lighthill, M. J. : Oscillating Airfoils at High Mach Number, Journal of 
the Aeronautical Sciences, Vol. 20, No. 6, June 1953. 

25. Jones, R. T. : The Unsteady Lift of a Wing of Finite Aspect Ratio, NACA 
Report 681, 1940. 

26. Moler, C. B. and G. W. Stewart: An Algorithm for Generalized Matrix 
Eigenvalue Problems, SIAM Journal of Numerical Analysis, Vol. 10, No. 2, 
April 1973. 

27. Ward, R. C. : An Extension of the QZ Algorithm for Solving the Generalized 
Matrix Eigenvalue Problem. NASA-TN-D7305, July 1973. 

28. Smith, B. T. , et al : Lecture Notes in Computer Science, No. 6, Matrix 
Eigensystem Routines - EISPACK Guide. Springer-Verlag, Berlin, 1976. 

29. Garbow, B. S., et al : Lecture Notes in Computer Science, No. 51, Matrix 
Eigensystem Routines - EISPACK Guide Extention. Springer-Verlag, Beilin, 
1977. 

30. Roark, R. J. and W. C. Young: Formulas for Stress and Strain. McGraw-Hill 
Book Co., New York, 1975, pp. 294. 

31. Gessow, A. and G. C. Myers: Aerodynamics of the Helicopter. Frederick 
Ungar Publishing Co., New York, 1967 

32. Kuczynsky, W. A.: Inflight Rotor Stability Monitor. NASA Symposium on 
Flutter Testing Techniques, NASA Flight Research Center, Edwards AFB, 
California, October 1975. 

33. Bielawa, R. L. : Dynamic Analysis of Multi-Degree-Of-Freedom Systems Using 
Phasing Matrices. American Helicopter Society/NASA Specialists Meeting 
on Rotorcraft Dynamics Proceedings, No. 4, Moffett Field, California, 

1974. 

34. Hall, W. E. : Application of Floquet Theory to the Analysis of Rotary Wing 
VTOL Stability. SUDDAR No. 400, Stanford University Center for Systems 
Research. (NASA Contract NAS-25143), 1973. 


248 



3. Recipient’s Catalog No. 



Aeroelastic Analysis for Propellers - Mathematical 
Formulations and Program User's Manual 

7. Authors) 

R. L. Bielawa, S. A. Johnson, R. M. Chi, and 

S . T . Gangwani 

9. Performing Organization Name and Address 

United Technologies Research Center 
East Hartford, Connecticut 06108 


5. Report Date 

December 1983 

6. Performing Organization Code 

8. Performing Organization Report No. 

UTRC83-6 

To. Work Unit No. 

11. Contract or Grant No. 

NAS3-22753 

13. Type of Report and Period Covered 


12. Sponsoring Agency Name and Address Contractor Report 

National Aeronautics and Space Administration u. Sponsoring Agency Code 

Washington, D.C. 20546 535-03-12 (E-1799) 

15. Supplementary Notes 

Final report. Project Manager, Oral Mehmed, Structures and Mechanical Technologies 
Division, NASA Lewis Research Center, Cleveland, Ohio 44135. 

16. Abstract 

Mathematical development is presented for a specialized propeller dedicated version 
of the United Technologies Corporation G400 Rotor Aeroelastic Analysis. This 
specialized analysis, G400PROP, simulates aeroelastic characteristics particular 
to propellers such as structural sweep, aerodynamic sweep and high subsonic 
unsteady airloads (both stalled and unstalled) . Detailed formulations are 
presented for these expanded propeller related methodologies. Results are 
presented of limited application of the analysis to realistic blade configurations 
and operating conditions which include stable and unstable stall flutter test 
conditions. Sections are included for enhanced program user efficiency and 
expanded utilization. This material includes (1) a detailed description of the 
structuring of the G400PROP FORTRAN coding, (2) a detailed description of the 
required input data, (3) a detailed description of the output results, and (4) 
general information to facilitate operation and improve efficiency. 


i 


17. Key Words (Suggested by Authors)) 

Propellers 

Rotor aeroelasticity 
Structural sweep 
Unsteady airloads 
Stall flutter 


18. Distribution Statement 
Unclassified 

STAR Category 

- Unlimited 
39 


19. Security Classif. (of this report) 

20. Security Classif. (of this page) 

21. No. of pages 

22. Price* 

Unclassified 

Unclassified 

253 

A12 


*For sale by the National Technical Information Service, Springfield, Virginia 22161 


NASA-Langley, 1983 









